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1.  Random  Number  Generators 


The  random  number  generators  (RNGs)  evaluated  in  this  report  fall  into  a  few  main  categories. 
The  RNG  algorithms  are  listed  in  appendix  A.  Initialization  (seeding)  algorithms  are  not 
presented. 

The  first  category  includes  variations  on  linear  congruential  generator  (LCG)  algorithms. 

Library  functions,  denoted  lib .  rand,  lib .  rand_r,  and  lib .  random,  are  generally  not  portable. 
These  are  executed  from  library  calls,  and  the  implementations  depend  on  platform.  The  extent 
of  standardization  (1 )  is  only  that:  “The  rand  function  computes  a  sequence  of  pseudo-random 
integers  in  the  range  0  to  RAND  JV1AX”  and  “The  value  of  the  RAND_MAX  macro  shall  be  at 
least  32767.” 

In  an  effort  to  be  specific,  cyg .  rand  is  from  Cygwin  C  library  source  code,  glib .  rand_r  is  from 
GNU  C  library  source  code,  and  Cstd .  randl  is  from  the  essential  reference  (2).  On  some 
platforms,  these  may  coincide  with  library  functions,  although  this  is  not  required  by  the  standard. 

The  MUVES-S2  RNG  rnrand  is  a  variant  of  Cstd.  randl. 

The  standard  Java  RNG  java32  is  from  Java  library  source. 

Multiplicative  linear  congruential  generator  (MLCG)  algorithms  include  the  Park- Miller  minimal 
standard  RNG  with  Bays-Durham  shuffle  (3)  pm3bds,  which  is  the  original  RNG  implementation 
in  MUVES  3,  and  cyg .  rand_r,  taken  from  Cygwin  C  library  source  code. 

The  next  two  groups  are  based  on  shift  register  (SR)  implementation. 

RNGs  based  on  the  linear  feedback  shift  register  (LFSR)  include  Tausworthe  ( 4 ,  5)  RNGs,  the 
32-bit  taus©88  and  tausll3,  and  the  64-bit  taus258  (tested  as  taus258x  and  taus258x32). 

Variations  on  generalized  feedback  shift  register  (GFSR)  algorithms  include  the  Ziff  ( 6 )  RNG 
gfsr,  a  four-tap  GFSR  algorithm,  the  Mersenne  Twister  (7,  8)  32-bit  version  mt  and  64-bit 
version  mt64  (tested  as  mtx  and  mtx32),  and  the  well-equidistributed  long-period  linear  (9)  RNG 
welll©24a. 

RNGs  due  to  Marsaglia  are  multiply  with  carry  (MWC)  variations  cmwc4©96,  mwc,  mwcx, 
exclusive-or  (XOR)  Shift  (10)  xorl28,  and  keep  it  simple  stupid  kiss,  a  combination  of  FFSR, 
FCG,  and  MWC  elements. 
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2.  Testing  and  Selection 


Major  test  suites  for  RNGs  are  TestUOl  (if),  the  NIST  suite  (72),  and  the  Diehard  suites  DH1 
and  DH2  (13).  The  Diehard  suites  implement  many  classical  RNG  tests,  the  NIST  suite  is  geared 
to  cryptographic  security,  and  TestUOl  is  extremely  diverse,  implementing  classical  tests, 
cryptographic  tests,  new  tests  proposed  in  the  literature,  and  original  tests. 

All  RNGs  were  tested  with  the  complete  DH1  and  DH2  suites  and  NIST  suite.  All  RNGs  except 
for  some  of  the  library  functions  and  algorithms  were  subjected  to  the  major  test  batteries  from 
TestUOl. 

Since  all  suites  use  32-bit  integers,  the  64-bit  algorithms  were  used  to  create  32-bit  output  either 
by  taking  the  upper  32  bits  in  the  case  of  mtx  and  taus258x  or  by  alternating  upper  and  lower 
32-bit  segments  in  the  case  of  mtx32  and  taus2 58x32. 

DH1  implements  24  tests  and  DH2  implements  26  tests.  Each  test  produces  a  single  p-value  p. 
We  regard  that  test  as  failed  if  p  <  0.005,  indicating  departure  from  uniformity. 

The  NIST  suite  implements  189  tests  and  reports  p’s.  A  test  is  considered  failed  if  p  <  0.0001  or 
p  >  0.9999. 

The  TestUOlOl  suite  implements  773  tests  and  reports  p’s.  A  test  is  considered  failed  if 
p  <  0.001  or p  >  0.999. 

The  p-value  selection  criteria  for  the  various  test  suites  were  chosen  to  produce  a  few  failures  in 
the  best  cases.  Setting  the  criteria  too  low  (closer  to  zero)  would  exhibit  no  failures,  and  setting 
the  criteria  too  high  would  fail  everything. 

Table  1  tallies  the  failure  counts  for  each  RNG  by  test  suite.  Detailed  test  results  are  available  on 
the  local  intranet  (14)  or  by  request  from  the  author. 

SR  (LFSR  and  GFSR)  RNGs  are  expected  to  fail  Finear  Complexity  and  Fempel-Ziv 
Compression  tests.  This  causes  concern  for  cryptographic  use  but  is  not  relevant  for  our 
applications.  In  fact,  the  features  of  these  algorithms  exploited  in  section  3  imply  such  failures. 

Additionally,  FFSR  (e.g.,  Tausworthe)  RNGs  are  known  to  fail  Matrix  Rank  tests.  Here,  each 
Tausworthe  RNG  fails  eight  TestUOl  Matrix  Rank  tests.  This  may  cause  problems  when 
simulating  large  binary  matrices.  Curiously,  well  1024a,  which  is  supposed  to  avoid  some  of  the 
problems  of  the  Mersenne  Twisters,  fails  four  TestUOl  Matrix  Rank  tests.  The  Mersenne 
Twisters  pass  all  TestUOl  Matrix  Rank  tests. 

The  NIST  suite  also  tests  Matrix  Rank.  Marginal  failure  for  mt32x  is  the  only  NIST  Matrix  Rank 
test  failure  among  the  FFSR  and  GFSR  RNGs. 
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Table  1.  RNG  test  suite  failure  counts. 


RNG 

Test  Suite 

Category 

Name 

(size) 

DH1 

DH2 

NIST 

TestUOl 

TestU01a 

LCG 

cyg . rand 

(32) 

1 

1 

7 

17 

16 

glib.rand_r 

(31) 

16 

20 

147 

- 

- 

lib . rand 

(31) 

13 

15 

148 

- 

- 

lib . random 

(31) 

13 

15 

149 

- 

- 

lib .  rand_r 

(31) 

15 

18 

146 

- 

- 

Cstd.randl 

(31) 

18 

23 

148 

- 

- 

rnrand 

(15) 

17 

22 

187 

404 

375 

java32 

(32) 

4 

6 

9 

111 

109 

MLCG 

pm3bds 

(31) 

5 

10 

149 

157 

154 

cyg .  rand_r 

(31) 

15 

18 

146 

- 

- 

LFSR 

taus088 

(32) 

0 

0 

4 

15 

1 

tausll3 

(32) 

0 

0 

5 

16 

2 

taus258x 

(64) 

0 

0 

6 

12 

0 

taus258x32 

(64) 

0 

1 

6 

13 

1 

GFSR 

gfsr 

(32) 

0 

0 

6 

4 

0 

mt 

(32) 

1 

1 

7 

6 

2 

mtx 

(64) 

1 

0 

2 

4 

0 

mtx32 

(64) 

1 

0 

9 

5 

1 

welll024a 

(32) 

0 

0 

11 

9 

1 

MWC 

cmwc4Q96 

(32) 

0 

1 

3 

1 

1 

mwc 

(32) 

0 

1 

8 

2 

2 

mwcx 

(32) 

0 

0 

8 

2 

2 

XOR  Shift 

xorl28 

(32) 

0 

0 

9 

30 

16 

Combined 

kiss 

(32) 

0 

0 

5 

2 

2 

a  Linear  Complexity,  Lempel-Ziv  Compression,  and  Matrix  Rank  test  failures  excluded. 


The  Mersenne  Twister  is  widely  accepted  by  the  community;  in  fact,  mt  is  the  core  RNG 
implementation  in  the  Python  language  and  is  a  default  RNG  in  Mathematica,  PHP,  and  GLib. 

The  Tausworthe  RNGs  are  efficient  and  can  be  extremely  fast. 

Based  on  these  considerations  and  the  test  results,  the  author  recommends  use  of  the  32-bit 
Mersenne  Twister  mt  as  a  primary  RNG  and  the  four-SR  Tausworthe  tausll3  as  secondary.  If 
64-bit  algorithms  are  required,  the  64-bit  Mersenne  Twister  mt64  and  five-SR  Tausworthe 
taus258  algorithms  should  be  implemented. 

These  RNGs  share  a  common  structure  that  facilitates  the  implementation  of  statistically 
independent  RNGs  as  shown  in  the  following  section. 
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3.  F2-Linear  RNGs 


3.1  Motivation 

A  wide  class  of  fast  long-period  RNGs  are  based  on  linear  recurrence  modulo  2.  This  class 
includes  the  LFSR  Tausworthe  and  the  GFSR  Mersenne  Twister  family  RNGs.  The  RNGs 
considered  here  have  approximate  periods  p  ranging  from  2 88  3  ■  1026  to  219937  ~4  -  lO6001. 

A  necessary  feature  of  RNG  implementation  is  the  availability  of  statistically  independent  RNGs. 
The  creation  of  algorithmically  independent  RNGs  is  trivial.  One  simply  allocates  separate  state 
information  for  each  instance  of  the  RNG  and  assures  each  RNG  modifies  only  its  own  state. 

However,  this  does  not  assure  statistical  independence.  Consider  that  any  RNG  generates  a  single 
output  sequence,  and  selecting  the  initial  state  of  the  RNG  (called  seeding)  chooses  a  particular 
starting  place  in  the  sequence.  Two  instances  are  considered  to  be  independent  if  their  sequences 
do  not  overlap  during  execution  of  the  program  using  the  RNGs. 

To  obtain  two  maximally  independent  instances  of  an  RNG,  one  could  choose  two  starting  points 
separated  by  p/2.  One  way  to  do  this  is  to  seed  the  first  RNG,  start  the  second  RNG  at  the  same 
point,  and  execute  the  second  instance  p / 2  times  to  obtain  its  start  state.  The  simplest 
Tausworthe  RNG  has  period  p  ~  288.  One  would  initialize  the  second  instance  by  executing  it 
287  times.  If  the  RNG  executes  1  billion  or  109  ~  230  times/second,  this  would  take  about  257  s. 
Since  there  are  about  225  s  in  a  year,  the  RNGs  would  be  ready  to  use  in  about  232  ~  4  •  109  years. 
That’s  4  billion  years.  The  other  Tausworthe  RNGs  to  consider  have  periods  2 113  and  2258. 

The  Mersenne  Twister  has  p  ~  2 19937  —  1.  Initializing  it  in  this  way  would  take  >  2  ■  105991  years. 
This  is  not  practical. 

Fortunately,  the  technology  shown  here  can  get  this  down  to  microseconds  for  the  Tausworthe 
RNGs  and  a  few  milliseconds  for  the  Mersenne  Twister.  The  methods  are  due  to  L‘Ecuyer  et 
al.  (15,  16). 

3.2  F2-Linear  RNG  Algorithms 

All  the  RNGs  mentioned  previously  work  in  essentially  the  same  way.  The  theory  is  based  on 
arithmetic  in  the  finite  field  with  two  elements,  F2  —  {0, 1}.  The  binary  operations  of  addition 
and  multiplication  in  F2  are  characterized  by 

0  +  0  =  1  +  1  =  0,  and  0+  1  =  1  +  0  =  1  (1) 

and  0x0=lx0  =  0xl  =  0,  and  1x1  =  1.  (2) 
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Since  1  is  its  own  additive  inverse,  a  —  —a  in  F2,  and  subtraction  is  the  same  as  addition. 
Bitwise  “exclusive  or”,  denoted  ',  is  addition,  so  a+b  =  a~b.  Bitwise  “and”,  denoted  &, 
is  multiplication,  so  a  x  b  =  a&b. 

State  information  for  the  RNG  is  held  in  x,  an  element  of  the  w -dimensional  coordinate  vector 
space  F2  .  In  simple  cases,  x  is  an  unsigned  integer,  and  w  =  32.  The  RNG  is  implemented  by 
recursive  linear  transformation  on  the  state 


xn+\  =  Bxn ,  where  n  =  0, 1, 2, . . .  (3) 

to  give  the  RNG  sequence  x0,xi,x2, ....  Multiple  steps  are  realized  through  repeated  application 
of  the  transformation,  symbolized  by  powers  of  B  as  xn  —  Bnx 0  or,  in  general,  xn+i  =  B"xt.  The 
identity  transformation  is  B°  —  I. 

All  the  polynomials  we  consider  are  elements  of  F2[z],  the  ring  of  polynomials  in  indeterminate  z 
with  coefficients  in  F2.  The  transformation  B  has  a  characteristic  polynomial  (CP)  of  degree 

k  ^  w, 

C(z)  =  det(B-zI)  =  zk  +  ck_izk  1  +  •  •  • +  ciz +  co,  (4) 

which  satisfies  C(B)  —  0  according  to  the  Cayley-Hamilton  theorem  (any  linear  transformation  is 
a  zero  of  its  own  CP).  Therefore  Bk  can  be  expressed  as  a  combination  of  lower-power  terms, 

Bk  —  Ck-\Bk  1  +  •  •  •  -\-C\B  -\-cqI.  (5) 

To  compute  some  future  state  x(i  directly  from  the  initial  state  jto, 

xd  =  Bdx  o,  (6) 

apply  the  division  algorithm  to  the  polynomial  zd  to  obtain 

zd  =  Q(z)-C{z)  +  J{z),  (7) 

where  /  =  0  or  deg  J  <  deg  C.  We  write 

J{z)  —  zdmod C(z)  —  jk-\Zk  1  +jk-2Zk  2H - f  jiz2  +jiZ+jo-  (8) 

Since  C(B)  —  0,  we  have  Bd  —  J(B),  and  the  transformation  to  state  d  is  an  TVlinear 
combination  of  the  transformations  Bk~l,  ■■■  ,B°  =  I,  with  /,  e  {0,  1}, 

Bd=jk_iBk-1  +  -+j0B°,  (9) 

and  state  xd  is  the  same  C2 -linear  combination  of  the  k  initial  states  x-t  =  B'x o, 

xd  —  jk-\xk-l  H - bj o^o.  (10) 
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where  arithmetic  in  the  coordinate  vector  space  F T  yields  the  final  result,  the  jump  state  xj. 

The  CP  C  is  only  calculated  once  per  RNG.  The  jump  polynomial  (JP)  J  must  be  calculated  for 
each  d.  But  this  can  be  stored  and  applied  to  any  RNG  state  to  obtain  the  state  d  iterations 
forward.  The  Tausworthe  generators  use  three,  four,  or  five  parallel  32-bit  or  64-bit  words  and 
thus  have  w  ^  64  and  require  fewer  than  64  states  to  perform  the  calculation  with  no  more  than 
five  words  of  state  information  (this  is  accomplished  in  microseconds).  The  Mersenne  Twisters 
use  624  32-bit  words  or  312  64-bit  words  for  19,968  bits  of  state  information.  These  require 
fewer  than  21,000  states  for  calculation  of  any  future  state  and  executes  in  milliseconds. 

3.3  An  Example 

Consider  the  first  subgenerator  of  taus088.  The  algorithm  is  denoted  gen,  and  calling  gen() 
updates  x  to  B.x.  In  generic  C/C++  code,  the  RNG  algorithm  follows. 

unsigned  x  =  ®xl2B40fQQ; 
unsigned  gen  (  void  )  { 

x  =  (  (  x  &  Sxfffffffe  )  «  12  )  "  (  (  (  x  «  13  )  "  x  )  »  19  ); 
return  x; 

} 

The  initial  (seed)  state  of  x  is  given.  Executing  gen()  advances  x  to  the  next  state  and  returns 
that  value  to  be  used  as  the  random  number.  Details:  x  is  a  32-bit  unsigned  integer,  &  is  bitwise 
and,  «  is  shift  left,  »  is  shift  right,  and  '  is  bitwise  xor.  Note  that  &  is  multiplication  and  "  is 
addition  (mod  2). 

The  matrix  of  the  transformation  is  never  used  in  computation. 

C  and  J  are  stored  in  32-bit  unsigned  integers  where  each  bit  is  a  coefficient,  and  the  constant 
term  is  the  least  significant  bit.  So  c,  =  C  »  i  &  1  and  y,  =  1  »  i  &  1. 

The  CP,  given  by  the  Berlekamp-Massey  algorithm  (17,  18),  is 

C(z)=z31+z25+z19+z13  +  l,  (11) 

and  its  coefficients  are  0x82082001  =  1000  0010  0000  1000  0010  0000  0000  0001 2. 

For  a  jump  displacement  of  d  —  220,  we  have  J  =  zn  mod  C  =  z2'"  mod  C  =  z1048576  mod  C.  The 
JP  is 


J(Z)  =  z30  +  z27  +  z26  +  z25  +z24  +  z23  +z21  +  z2U  +z19  +z 


,26 


,25 


,23 


,20 


T9 


,18 


+  z14  +  z12+z9+z8+z5, 


(12) 
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and  its  coefficients  are  J  —  0.v4/£c5320  =  0100  1111  1011  1100  0101  001 1  0010  00002,  so 
the  jump  state  x,i  =  *1048576  is  the  same  F2 -linear  combination  of  base  states  in  F™, 


* 1048576  —  *30  +*27  +*26  +*25  +*24  +*23  +*21  +*20  +*19  +*18 

+  *14  +*12  +  *9  +  *8  +  *5  •  (13) 

The  basic  jump  state  calculation  algorithm  follows.  To  calculate  x  —  x,j  in  state  d,  start  with  x  in 
an  initial  state  x  =  jc0,  and  collect  the  32  states  *0  =  5  [0],  . . .  ,  jc3i  =  s  [3 1]  by 


for  (  i=8  ;  i<31  ;  i++)  { 
s [i]  =  x; 
genO  ; 

} 


Then  calculate  x  =  xj  by 


bitset<32>  J  =  0x4fbc5320; 
x  =  0; 

for  (  i=®  ;  i<32  ;  i++) 
if  C  JCi]  ) 
x  ~=  s [i]  ; 


3.4  Computation  Details 

A  complete  C++  implementation  of  this  system  for  five  RNGs  is  presented  in  appendix  B.  The 
RNGs  are  demo  (the  first  SR  of  t®88  as  in  the  example  of  section  3.3),  tll3  (the  32-bit  four-SR 
Tausworthe  RNG),  t258  (the  64-bit  five-SR  Tausworthe  RNG),  mt32  (the  32-bit  Mersenne 
Twister),  and  mt64  (the  64-bit  Mersenne  Twister).  Extract  the  files  of  appendix  B  to  a  common 
directory  and  make.  The  source  code  is  available  on  the  local  intranet  (14)  or  by  request  from  the 
author. 

Each  resulting  executable  performs  CP  and  JP  calculations  for  one  of  the  RNGs.  (The  CP  results 
have  been  inserted  into  the  code  for  JP  calculations.)  The  driver  poly_gen .  cpp  contains  usage 
instructions. 

Computational  details  are  presented  in  sections  3. 4. 1-3 .4. 4. 
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3.4.1  Recurrence  Relation  for  F2 -Linear  RNGs 

Efficient  calculation  of  the  CP  relies  on  the  following  recurrence  property.  Since  C(B)  —  0  and 
Bk  —  cu-\Bk~l  H - \- c \B  -\- cqI ,  the  RNG  is  also  defined  by 

Xk  —  Bkx0  —  (cfc-\Bk  1  +  •  ■  •  -\-C\B  +co/)^o 

—  Ck-\Bk  lXo  +  •  •  •  +  C  i  Bxq  +  CqIXq 

—  Ck-iXk-i  4 - \-c\X\  +co^o>  (14) 

leading  to  the  recurrence  relation  for  all  i  ^  0, 

Xk+i  =  Ck-lXk-l+i  H - h  CqXi  ,  (15) 

or  equivalently  for  all  n  ^  k, 

X n  Ck—lXfi—i  ~h  ‘  ‘  '  ~h  C()Xn— k  •  (16) 

Obviously,  each  bit  of  x  obeys  the  same  relation. 

3.4.2  The  Characteristic  Polynomial 

Apply  the  Berlekamp-Massey  algorithm  (17,  18)  to  a  single  bit  stream  from  the  RNG.  This 
computes  the  minimal  polynomial  of  the  linearly  recurrent  sequence  in  equation  16  and  the  RNG 
(equation  3).  If  the  minimal  polynomial  is  primitive,  then  it  is  the  CP. 

The  stream  must  have  length  exceeding  twice  the  degree  of  the  desired  CP.  Any  bit  position  will 
do  since  they  all  satisfy  the  same  recurrence.  The  algorithm  operates  on  a  polynomial  whose 
coefficients  are  given  by  the  bit  stream.  The  implementation  in  appendix  B  obtains  the  degree 
19937  CPs  of  the  32-bit  and  64-bit  Mersenne  Twisters  in  less  than  4  s  each  on  a  Windows 
Pentium  desktop  computer  and  less  than  2.5  s  each  on  a  Linux  Xeon  server. 

They  have  the  same  degree,  but  the  32-bit  Mersenne  Twister  characteristic  polynomial 

C(z)=z19937  +z19314  +  z19087  +z18860  +  ---+z1416  +  z1189  +  l  (17) 

has  135  nonzero  coefficients,  and  the  64-bit  Mersenne  Twister  CP 

C(z)  =  z 19937  +  z19626  +  z 19470  +  z 19314  +  •••  +  z468  +  z312  +  1  (18) 

has  285  nonzero  coefficients. 

Tausworthe  RNG  CPs  are  calculated  in  less  than  0.2  s  on  the  Pentium  and  less  than  0.04  s  on  the 
Xeon.  Verification  of  polynomial  primitivity  was  not  performed. 


8 


3.4.3  The  Jump  Polynomial 


The  JP  for  jump  size  d  is 


J  (z)  =zd  mod  C(z),  (19) 

where  C  is  the  CP.  A  clever  algorithm  of  Knuth  (19)  allows  calculation  of  J  for  arbitrarily  large 
d  in  log  J  time  using  polynomials  of  degree  C  2  deg  C.  Obtain  f(z)  —  zd  from  the  sequence 

fo(x)  =z,  ...  ,  fk(z)  =zn,  (20) 


where 


k 

rf=5>-2‘-‘ 

(=0 


(21) 


and  the  binary  representation  of  d  is  given  by  the  sequence  (bo,  b  i,  . . . ,  bk  )  with  b{)  =  1  and 
bj  G  (0,  1}  for  i  >  0.  The  sequence  /,•  is  defined  recursively  by  fo(z)  =  z  and  for  1  ^  i  k. 


m 


Z  ■  fi-i(z)2 ,  bj  =  1 
fi~i(z)2,  bj  —  0. 


(22) 


For  example,  with  d 


876  = 


z ,  z 


11  0110  1 IOO2,  we  have  k  =  9  and  the  sequence  is 

_6  „13  21  54  109  „219  ^438  ^876 

Z  >  Z  y  Z  y  Z  y  Z  y  Z  y  Z  y  Z 


(23) 


At  each  step,  the  previous  value  is  squared  if  b  —  0  or  squared  and  multiplied  by  z  if  b  =  1 .  Proof 
is  by  induction  on  k.  With  k  =  1,  we  have  either  d  —  2=  IO2  and  the  sequence  is  z,  z2  or 
d  —  3  =  1 12  and  the  sequence  is  z,  z3.  Suppose  the  proposition  holds  for  length  k,  and  consider  d 
with  a  binary  sequence  one  larger: 


d  =  (b0,bi,...,bk,bk+i),  (24) 

where  dQ  =  (b0,bu...,bk).  (25) 

Either 

d=  (bo,bi,...,bk,  0)  =  2-d0  (26) 

or  d=  (bQ,bi,...,bk,  1)  =  2-d0+  1.  (27) 


In  the  first  case  zd  —  ( zdo )2,  and  in  the  second  case  zd  —  z-  ( zd° )2.  Since  zdo  is  obtained  by  the 
induction  hypothesis,  the  assertion  is  established. 
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To  compute  the  JP  /  =  zn  mod  C,  apply  the  algorithm  modulo  C(z)  as  follows: 


fi(z)  modC(z) 


z  ■  fi-i(z)2  modC(z) ,  bi=  1 
/•_i(z)2modC(z),  b\  —  0 


(28) 


3.4.4  Polynomial  Arithmetic 

Jump  polynomial  calculation  requires  evaluation  of  z  •  a(z)  and  ci(z  ')2  in  the  quotient  ring 
Q  =  F2[z]/C(z).  Let  «  =  deg(C).  Any  a(z)  =  a  a1  G  2  has  deg(a)  <  n.  If  a„_i  =  0,  then 

z-a(z)  still  lies  in  Q.  If  a„_i  =  1,  the  remainder  from  long  division  of  z  •  a(z)  by  C(z)  is  the 
result  z  -a(z)  mod  C{z)  —z-a(z)  —  C(z). 

Any  a  G  Q  can  be  represented  in  a  shift  register  a  of  length  n  with  the  constant  term  a{)  on  the 
right  and  an_i,  the  coefficient  of  zn  1 ,  on  the  left.  Multiplication  by  z  is  the  shift-left  operation 
«.  Subtraction  is  addition,  implemented  by  “exclusive  or”,  denoted 

This  algorithm  replaces  a  =  a(z)  by  z-a(z)  mod  C(z)  for  either  value  of  a„_ 

x  =  a[n-l] ; 
a  «=  1 ; 

if  (  x  )  a  ~=  c; 

A  similar  algorithm  replaces  a  =  a(z)  by  a(z)2 mod C(z). 
t.resetC) ; 

for  (  int  j=n-l;  j>=Q;  j--  )  { 
x  =  t [n-1] ; 
t  «=  1; 

if  (  x  )  t  "=  c; 
if  (  a[j]  )  t  '=  a; 

} 

a  =  t; 

Here,  t  is  another  shift  register  of  length  n ,  and  t .  reset  ()  sets  all  bits  of  t  to  0.  This  algorithm 
is  an  explicit  implementation  of  the  product  a(z)  ■ b(z )  as 

n—  1 

a  b  —  '^2(aj-zj -b)  =  a0b  +  z{a\b  +  z(a2b  +  z( - \-z(an_2b  +z(an^lb))  •••))).  (29) 

7=0 

3.5  Implementation 

Some  care  must  be  exercised  in  the  sequence  of  seeding,  RNG  algorithm  execution  for 
incrementing  the  state,  and  state  information  collection  for  jump  calculation.  For  the  purpose  of 
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jump  calculations,  the  main  concern  is  the  state  of  the  RNG  and  not  the  random  integer  output 
(which  is  a  side  effect  of  state  transition). 

Seed  states  are  not  useful  for  jump  calculations,  and  the  calculations  can  fail  if  a  seed  state  is 
included  in  equation  10.  The  remedy  for  this  is  to  ensure  that  all  states  used  in  jump  calculations 
are  given  by  RNG  state  transition.  In  the  following  implementations,  a  single  call  to  the  RNG 
algorithm  itself  is  included  in  any  seeding  procedure  so  that  the  RNG  is  in  a  valid  algorithm  state, 
called  the  zero  state  and  denoted  xq. 

A  collection  of  statistically  independent  FVlinear  RNGs  can  be  implemented  as  follows.  Seed 
the  first  RNG  arbitrarily  (and  execute  the  RNG  algorithm  once  as  the  previous  paragraph 
indicates).  Save  the  resulting  zero  state  as  a  reference  state.  Instantiate  subsequent  RNGs  by 
applying  the  jump  procedure  to  the  reference  state,  setting  the  state  of  the  new  RNG  to  the  jump 
state  and  saving  the  jump  state  as  the  new  reference  state. 

Calculation  of  the  CP  is  done  once  and  for  all  for  each  RNG.  The  CP  is  used  to  compute  the  JP, 
and  the  CP  does  not  appear  in  the  implementation.  A  single  jump  displacement  d  serves  as  a 
common  increment  for  the  system.  Therefore,  calculation  of  the  JP  need  be  done  once  for  the 
given  displacement  and  can  be  applied  to  any  state  configuration  to  calculate  a  jump  state. 

Using  the  zero  state  as  a  reference  state  for  jump  calculations,  the  state  sequence  is  xq,x\,x2,  . 
and  jump  to  state  xd  is  calculated  as 


xd  —  jk-iXk-1  H - \-  j  cAo>  (30) 

where  equation  8  gives  the  JP  coefficients  (For  small  d  with  J  —  zd,  we  have  /,■  =  1  if  i  —  d 
and  ji  —  0  otherwise,  giving  jump  state  xd  as  required.)  From  this  point,  the  state  sequence  is 
xd,Xd+i,xt i+2, ....  Using  xd  as  a  reference  state  gives  jump  state  X2d,  from  which  point  the  state 
sequence  is  x2d,  x2d+i,x2d+2, 

There  is  a  tradeoff  between  the  jump  increment  d  and  the  number  k  of  available  RNGs.  Each 
RNG  will  run  for  d  iterations  before  it  overlaps  the  next  RNG.  If  the  period  of  the  base  RNG  is  p, 
the  relation  is 

p=kd,  (31) 

so,  effectively,  k—p/d  independent  RNGs  of  “period”  length  d  are  available. 

A  complete  C++  implementation  of  this  system  for  five  RNGs  is  presented  in  appendix  C.  The 
RNGs  are  demo  (the  first  SR  of  t®88  as  in  the  example  of  section  3.3),  tll3  (the  32-bit  four-SR 
Tausworthe  RNG),  t258  (the  64-bit  five-SR  Tausworthe  RNG),  mt32  (the  32-bit  Mersenne 
Twister),  and  mt64  (the  64-bit  Mersenne  Twister).  Extract  the  files  of  Appendix  C  to  a  common 
directory  and  make.  The  source  code  is  available  on  the  local  intranet  (14)  or  by  request  from  the 
author. 

Each  executable  illustrates  seeding,  allocation,  and  jump  displacement  verification  for  the 
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indicated  RNG.  JP  values  for  d  —  220  are  in  the  implementation  code  along  with  some  realistic 
values  that  could  be  used  in  practice.  Verification  is  performed  by  allocating  two  RNGs, 
executing  the  first  one  d  times,  and  checking  that  they  are  in  the  same  state. 

For  each  RNG,  the  procedure  was  tested  with  various  powers  of  2  from  21  to  approximately  235 
and  powers  of  10  from  101  to  10 10  or  so,  along  with  other  arbitrary  values.  Verification  of  larger 
d  values  is  time-consuming.  Verification  for  practical  values  of  d  is  impossible. 

This  system  generates  only  (32-bit  or  64-bit)  random  integers.  Other  discrete  and  continuous 
distributions  can  be  implemented,  for  example,  by  the  methods  of  Saucier  (20). 

3.5.1  Demo 

This  minimal  example  illustrates  the  basic  structure  and  operation  of  the  class  without  the  added 
detail  required  for  the  other  RNGs.  Calculations  proceed  as  in  section  3.3. 

3.5.2  Tausworthe 

For  the  Tausworthe  RNG  tll3,  we  have p  ~  21 13.  The  choice  of  d  —  280  implies  k  =  233  ~  1010 
independent  RNGs.  Generating  at  a  rate  of  230  ~  109/s,  the  RNGs  will  overlap  in  250  s,  or 
225  ~  3.2  x  107  years,  since  there  are  225  s  in  a  year.  This  should  give  an  adequate  run  time  and 
number  of  RNGs.  The  JP  for  d  =  280  is  in  the  implementation  code  but  commented  out. 

This  RNG  (see  section  A. 3,  page  23)  is  composed  of  four  independent  32-bit  SRs  that  do  not 
interact.  They  are  only  combined  to  produce  the  return  value  random  integer  which  has  no  part  in 
state  calculation.  Each  SR  has  its  own  CP  (calculated  using  a  bit  stream  from  that  SR)  and  its 
own  JP  for  any  given  jump  displacement.  Each  of  the  four  JPs  needs  32  bits,  and  these  are 
allocated  in  the  implementation  file  in  section  C.4  as 

unsigned  int  Jump_P[4]  =  { 

0x®c382e31  ,  ®xlb®4®425  ,  ®x®b49a5®9  ,  ®x®173f6b® 

}; 


and  then  the  jump  calculation  is  applied  separately  to  each  SR  via 

for  (  i  =  ®  ;  i  <  32  ;  i++  )  { 
for  (  j  =  ®  ;  j  <  4  ;  j++  ) 
if  (  JP[j]  [i]  ) 

temp_state.z[j]  "=  s.z[j]; 
gen() ; 

} 

where  each  J P [  j ]  is  a  bitset<32>  version  of  Jump_P[j] .  Corresponding  code  for  the  64-bit 
Tausworthe  RNG  t258,  with  five  independent  64-bit  SRs,  is  in  section  C.5,  with  JPs  for  d  =  280 
and  d  —  2 128  included  and  commented  out. 
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3.5.3  Mersenne  Twister 


For  the  Mersenne  Twister  mt32,  p  ~  219937,  and  d  =  2 10000  implies  k  =  2  9937  RNGs.  The  JP  for 
d  =  2 10000  is  in  the  implementation  code  but  commented  out. 

State  information  for  mt32  consists  of  the  structure  in  C.6: 


static  const  unsigned  N  =  624; 
struct  state  {  //  state  information 

unsigned  n; 
unsigned  y; 
unsigned  z [N] ; 

}  ; 


where  1  ^  n  N  and  n  —  1  is  a  static  index  into  the  state  array  z  with  y  =  z[n  —  1].  Explicit 
inclusion  of  y  in  the  state  structure  facilitates  state  collection  and  jump  calculation.  Upon  calling 
the  mg,  if  n  —  N,  the  array  z  is  recalculated  and  n  is  set  to  zero;  y  is  replaced  with  z[n\;  n  is 
incremented;  and  a  tempering  transformation  of  y  is  returned  as  RNG  output. 

Upon  seeding,  z  is  populated  in  some  fashion  (but  not  by  the  RNG  algorithm  itself),  and  n  is  set 
to  N.  The  first  RNG  execution  (call  0,  included  in  the  seeding  procedure)  then  computes  z  by  the 
RNG  algorithm,  sets  the  zero  state  y  =  z  [0] ,  and  sets  n  —  1.  The  next  N  —  1  calls  (call  i  with 
1  ^  i  ^  N  <- 1)  set  y  =  z\i\  and  n  =  i  +  1,  finally  exhausting  z  and  setting  n  —  N.  Then  the 
process  repeats,  recalculating  z  and  setting  y  =  z[0]  on  calls  Nk  for  integer  k. 

Two  RNG  instances  rO  and  rl  are  in  the  same  state  if  they  produce  the  same  RNG  output  stream 
or,  equivalently,  if  they  produce  the  same  y  sequence.  A  necessary  and  sufficient  condition  for 
state  equality  is  that  the  two  instances  generate  the  same  y  sequence  of  length  N.  Of  course,  this 
is  so  if  rO.n  =  rl.n  and  rO.y  =  r  l.y  and  rO.z  —  rl.z,  but  all  that  is  required  is  that  rO.y  =  rl.y 
and  the  next  N  —  1  values  of  y  agree.  Consider  that  the  period  p  —  2 19937  —  1  is  prime  and  that 
call  p  sets  the  same  state  y  as  call  0,  but  p  is  not  a  multiple  of  N  so  the  values  of  n  must  differ. 
Effectively,  z  is  a  moving  window  into  the  state  stream  y,  and  two  RNG  instances  can  generate 
the  same  stream  even  if  their  arrays  z  are  out  of  phase.  Any  state  structure  value  is  equivalent  to  a 
normalized  state  structure  with  n  =  1  and  y  =  z [0] .  The  array  z  is  an  artifact  of  computation. 

Usable  state  information  consists  of  the  state  sequence  y.  The  CP  is  calculated  from  a  single  bit 
position  ofay  sequence  of  length  >  2  ■  19937,  and  for  fixed  jump  displacement,  JP  is  calculated 
as  usual.  The  sequence  of  states  y  is  collected  as  a  basis  for  jump  calculation,  and  new  values  of 
y  are  computed  and  inserted  into  z  to  define  the  normalized  jump  state  structure.  Jump  state 
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calculation  requires  computing  N  jump  states  by 


yd  —  joyo+jiyi  H - \-jmym 

yd+i  —  joy  i  +  jiyi  H - \-jmyi+m 


yd+i  —  jayt + j  iy/+i  H - \-jmyi+m 


yd+N-l  —  joyN-l  +jiyN  ^ - \-  jmyN-l+m,  (32) 

where  m  —  deg (7)  <  19937  and  populating  the  array  z  with  the  jump  states.  Briefly, 

m 

ti  =  ^  jk  ■ yi+k  for  i  =  0, . . . ,  N  -  1 .  (33) 

k=0 

Then  set  n  =  1  and  each  z[i]  —  th  Calculations  are  always  based  on  a  reference  state  with  n  =  1, 
which  is  the  zero  state  (following  seeding  and  initial  RNG  call)  or  another  jump  calculation  state. 
No  more  than  19936  +  N  —  1  =  20559  states  y  must  be  collected  for  the  computation. 

Corresponding  code  for  the  64-bit  Mersenne  Twister  RNG  mt64  is  inlcuded  in  section  C.7,  with 
the  JP  for  d  =  2 10000  included  and  commented  out. 
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Appendix  A.  Random  Number  Generator  (RNG)  Algorithms  by  Category 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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Intentionally  Lett  Blank. 
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A.l  Linear  Congruential 
Cygwin  rand 

//  This  multiplier  was  obtained  from  Knuth,  D.E.,  "The  Art  of 
//  Computer  Programming,"  Vol  2,  Seminumerical  Algorithms,  Third 
//  Edition,  Addison-Wesley ,  1998,  p.  106  (line  26)  &  p.  108 

unsigned  CYG_RAND(void)  { 

seed  =  seed  *  6364136223846793005LL  +  1; 
return  (unsigned) (  (  seed  »  32  )  &  Oxffffffff) ; 


GNU  rand  r 

unsigned  GLIB_RAND_R(void)  { 
unsigned  int  next  =  Q; 
unsigned  int  result; 

next  *=  1103515245; 
next  +=  12345; 
next  &=  0x7fffffff ; 

result  =  (unsigned  int)  (next  »  16)  &  0x7ff;  //  0  . .  2047  =  2^11-1 

next  *=  1103515245; 
next  +=  12345; 
next  &=  0x7fffffff ; 
result  «=  10; 

result  *=  (unsigned  int)  (next  »  16)  &  0x3ff;  //  0  . .  1023  =  2~10-1 

next  *=  1103515245; 
next  +=  12345; 
next  &=  0x7fffffff ; 
result  «=  10; 

result  "=  (unsigned  int)  (next  »  16)  &  0x3ff;  //  0  .  .  1023  =  2~10-1 
Q  =  next; 

return  (unsigned)  result  &  Oxffffffff; 

} 
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C  Standard  randl 


unsigned  RAND 1 (void)  { 

Q  *=  1103515245 ; 

Q  +=  12345; 

Q  &=  0x7fffffff ;  //  2'31-1 

return  (unsigned)  Q  &  ©xffffffff; 

} 


MUVES-S2  rnrand 

unsigned  RNRAND (void)  { 

Q  *=  1103515245 ; 

Q  +=  12345; 

return  (unsigned)  (  Q  »  16  )  &  0x00007fff; 

} 


Java  java32 

unsigned  int  java32(  void  )  { 
int  bits=32; 

x  =  (x  *  0x5DEECE66Dull  +  OxBull)  &  ((lull  «  48)  -  1); 
return  (unsigned)  (  (  x  »  (48  -  bits)  )  &  Qxffffffffu  ) ; 

} 
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A.2  Multiplicative  Linear  Congruential 
Cygwin  rand_r 

//  Pseudo-random  generator  based  on  Minimal  Standard  by 
//  Lewis,  Goodman,  and  Miller  in  1969. 

// 

//  I[j  +  1]  =  a*I [j]  (mod  m) 

// 

//  where  a  =  16887  and  m  =  2147483647 

// 

//  Using  Schrage’s  algorithm,  a*I[j]  (mod  m)  can  be  rewritten  as: 

// 

//  a*(I[j]  mod  q)  -  r*{I[j]/q>  if  >=  8 

//  a*(I[j]  mod  q)  -  r*{I[j]/q}  +  m  otherwise 

// 

//  where:  {}  denotes  integer  division 
//  q  =  {m/a}  =  127773 

//  r  =  m  (mod  a)  =  2836 

// 

//  note  that  the  seed  value  of  8  cannot  be  used  in  the  calculation  as 
//  it  results  in  8  itself 

unsigned  CYG_RAND_R(void)  { 
int  k; 

int  s  =  (int) (seed); 
if  (s  ==  8) 

s  =  8x12345987; 
k  =  s  /  127773; 

s  =  16887  -  (s  -  k  -  127773)  -  2836  *  k; 
if  (s  <  8) 

s  +=  2147483647; 

(seed)  =  (unsigned  int)s; 
return  (unsigned) (s  &  Sxffffffff) ; 

} 
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Park-Miller  pm3bds 

#define  NTAB  32 

#define  M  0x7fffffff  //  214748B647  (Mersenne  prime  2 ^ 3 1- 1) 

#define  A  OxlOff5  //  69621 

#define  Q  0x787d  //  30845 

#define  R  Ox5d5e  //  23902 

static  int  next; 
static  int  DIV; 
int  table[  NTAB  ]; 

int  seed; 

unsigned  seed2 ; 

unsigned  PM3BDS(void)  { 

int  k  =  seed  /  Q;  //  seed  =  (  A*seed  )  %  M 

seed  =  A  *  (  seed  -  k  *  Q  )  -  k  *  R;  //  without  overflow  by 

if  (  seed  <  0  )  seed  +=  M;  //  Schrage’s  method 

int  index  =  next  /  DIV;  //  Bays-Durham  shuffle 

next  =  table  [  index  ] ;  //  seed  used  for  next  time 

table  [  index  ]  =  seed;  //  replace  with  new  seed 

return  next  «  1  ; 

} 


//  seed  to  be  used  as  index  into  table 

//  shuffle  table  of  seeds 
//  current  random  number  seed 
//  seed  for  tausworthe  random  bit 
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A.3  Linear  Feedback  Shift  Register 
Tausworthe  taus888 

#define  cl  Oxfffffffe  //  4294967294 
#define  c2  0xfffffff8  //  4294967288 
#define  c3  OxfffffffO  //  4294967280 

//  these  must  be  32-bit  integers 
static  unsigned  si; 
static  unsigned  s2; 
static  unsigned  s3; 

unsigned  taus088  (  void  )  { 

si  =  (  (  si  &  cl  )  «  12  )  “  (  (  (  si  «  13  )  ~  si  )  »  19  ) 

s2  =  (  (  s2  &  c2  )  «  4  )  "  (  (  (  s2  «  2  )  ~  s2  )  »  25  ) 

s3  =  (  (  s3  &  c3  )  «  17  )  *  (  (  (  s3  «  3  )  ~  s3  )  »  11  ) 

return  (si  s2  s3)  &  Oxffffffff; 

} 


Tausworthe  tausll3 

#define  cl  4294967294U 
#define  c2  4294967288U 
#define  c3  4294967288U 
#define  c4  4294967168U 

//  these  must  be  32-bit  integers 
static  unsigned  si; 
static  unsigned  s2; 
static  unsigned  s3; 
static  unsigned  s4; 

unsigned  tausll3  (  void  )  { 

si  =  (  (  si  &  cl  )  «  18  )  *  (  (  (  si  «  6  )  ~  si  )  »  13  ) 

s2  =  (  (  s2  &  c2  )  «  2  )  ~  (  (  (  s2  «  2  )  ~  s2  )  »  27  ) 

s3  =  (  (  s3  &  c3  )  «  7  )  "  (  (  (  s3  «  13  )  "  s3  )  »  21  ) 

s4  =  (  (  s4  &  c4  )  «  13  )  "  (  (  (  s4  «  3  )  ~  s4  )  »  12  ) 

return  (si  s2  ~  s3  s4)  &  Oxffffffff; 

} 
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Tausworthe  taus258 


#define  cl  OxFFFFFFFFFFFFFFFEULL  // 
#define  c2  0XFFFFFFFFFFFFFE0OULL  // 
#define  cB  0XFFFFFFFFFFFFF000ULL  // 
#define  c4  0XFFFFFFFFFFFE0000ULL  // 
#define  c5  0XFFFFFFFFFF800000ULL  // 


184467440737095 5 16 14ULL 
184467440737095 5 1104ULL 
18446744073709547520ULL 
1844674407  3  70942  0544ULL 
18446744073701 163008ULL 


//  these  must  be  64-bit  integers 
static  unsigned  long  long  si; 
static  unsigned  long  long  s2; 
static  unsigned  long  long  s3; 
static  unsigned  long  long  s4; 
static  unsigned  long  long  s5; 


unsigned  long  long  taus258  (  void  )  { 
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return  (si  s2  ~  s3  *  s4  s5)  &  QxffffffffffffffffULL; 

} 
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A.4  Generalized  Feedback  Shift  Register 
Generalized  Feedback  Shift  Register  gfsr 


static  const  int 
static  const  int 
static  const  int 
static  const  int 
static  const  int 


A 

B 

C 

D 

M 


=  471; 

=  1586; 

=  6988; 

=  9689; 

=  16383; 


//  period  is  2~D-1 
//  2" 14-1 


static  const  unsigned  long  int  MSB  =  0x80000000ul ; 
static  const  unsigned  long  int  MASK  =  ©xfffffffful ; 


static  const  int  WORD_SIZE  =  32 ; 


unsigned  long  int  _seed; 

unsigned  long  int  _state [16384] ;  //  _state[M+l]; 

int  _index ; 


unsigned  int  _int32(  void  )  {  //  return  an  integer  on  [0, 2 ~ 32-1] 
_index  =  (  _index  +  1  )  &  M; 

return  _state[  _index  ]  =  _state[  (  _index  +(M+1-A))&M] 
_state[  (  _index  +(M+1-B))&M] 

_state[  (  _index  +(M+1-C))&M] 

_state[  (  _index  +(M+1-D))&M]; 

} 
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Mersenne  Twister  mt32 


//  Period  parameters 
#define  N  624 
#define  M  397 

#define  MATRIX_A  0x9908b0df  //  constant  vector  a 
#define  UPPER_MASK  0x80000000  //  most  significant  w-r  bits 
#define  L0WER_MASK  0x7fffffff  //  least  significant  r  bits 

static  unsigned  mt[N];  //  the  array  for  the  state  vector 

static  int  mti=N+l; 

//  generates  a  random  number  on  [0 , Oxffffffff] -interval 
unsigned  mt32(void)  { 
unsigned  y; 
int  kk; 

static  unsigned  magOl [2] ={0x0 ,  MATRIX_A} ; 

if  (mti  >=  N)  {  //  generate  N  words  at  one  time 

for  (kk=0;kk<N-M;kk++)  { 

y  =  (mt[kk]&UPPER_MASK) | (mt[kk+l]&LOWER_MASK) ; 
mt[kk]  =  mt[kk+M]  *  (y  »  1)  *  mag01[y  &  0x1]; 

} 

for  ( ; kk<N - 1 ; kk++)  { 

y  =  (mt[kk]&UPPER_MASK) |(mt[kk+l]&LOWER_MASK) ; 
mt[kk]  =  mt[kk+(M-N)]  ~  (y  »  1)  ~  mag01[y  &  0x1]; 

} 

y  =  (mt [N- 1] &UPPER_MASK) | (mt [0] &LOWER_MASK) ; 
mt[N-l]  =  mt[M-l]  "  Cy  »  1)  "  mag01[y  &  0x1]; 

mti  =  0; 

} 

y  =  mt [mti++] ; 

//  Tempering 
y  “=  (y  »  11); 
y  "=  (y  «  7)  &  0x9d2c5680; 

y  “=  (y  «  15)  &  0xefc60000; 
y  "=  (y  »  18) ; 

return  y; 

} 
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Mersenne  Twister  mt64 


//  Period  parameters 
#define  NN  312 
#define  MM  156 

#define  MATRIX_A  0xB5026F5AA96619E9ULL 

#define  UM  0XFFFFFFFF8000O000ULL  //  Most  significant  33  bits 
#define  LM  0x7FFFFFFFULL  //  Least  significant  31  bits 

static  unsigned  long  long  mt[NN];  //  state  vector  array 
static  int  mti=NN+l; 

//  generates  a  random  number  on  [0,  2 "64-1] -interval 
unsigned  long  long  mt64(void)  { 
int  i ; 

unsigned  long  long  x; 

static  unsigned  long  long  magOl [2]={0ULL ,  MATRIX_A}; 

if  (mti  >=  NN)  {  //  generate  NN  words  at  one  time 

for  (i=0; i<NN-MM; i++)  { 
x  =  (mt[i]&UM) | (mt[i+l]&LM) ; 

mt[i]  =  mt[i+MM]  "  (x»l)  "  magOl  [(int)  (x&lULL)] ; 

} 

for  (;i<NN-l;i++)  { 
x  =  (mt[i]&UM) | (mt[i+l]&LM) ; 

mt[i]  =  mt[i+(MM-NN)]  "  (x»l)  "  mag01[(int) (x&lULL)] ; 

} 

x  =  (mt[NN-l]&UM) | (mt[0]&LM) ; 

mt[NN-l]  =  mt[MM-l]  "  (x»l)  "  magOl [(int) (x&lULL)] ; 


mti  =  0; 

} 

x  =  mt [mti++] ; 


//  Tempering 


X 

X 

II 

A 

A 

29) 

& 

0X5555555555555555ULL; 

X 

X 

II 

A 

A 
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& 

0x7 1D67FFFEDA60000ULL ; 

X 

X 
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II 
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A 
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0XFFF7EEEO00O00000ULL ; 
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return  x; 
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Well-Equidistributed  Long-Period  Linear  well  1824a 


#define  W  32 
#define  R  32 
#define  Ml  3 
#define  M2  24 
#define  M3  10 

#define  MATOPOS(t,v)  (v*(v»t)) 

#define  MAT0NEG(t,v)  (v* (v«(-(t)))) 

#define  Identity(v)  (v) 

#define  VO  STATE [state_i  ] 

#define  VM1  STATE [(state_i+Ml)  &  0x000000 lfU] 

#define  VM2  STATE [(state_i+M2)  &  0x000000 lfU] 

#define  VM3  STATE [(state_i+M3)  &  0x000000 lfU] 

#define  VRml  STATE [(state_i+31)  &  0x000000 lfU] 

#define  newVO  STATE [(state_i+31)  &  0x000000 lfU] 

#define  newVl  STATE [state_i  ] 

static  unsigned  int  state_i  =  0; 
static  unsigned  int  STATE [R] ; 
static  unsigned  int  zO,  zl,  z2; 

unsigned  WELLRNG1024a  (void)  { 
zO  =  VRml; 

zl  =  Identity(VO)  "  MATOPOS  (8,  VM1) ; 

z2  =  MATONEG  (-19,  VM2)  ~  MAT0NEG(- 14 , VM3) ; 
newVl  =  zl  z2; 

newVO  =  MATONEG  (-11, zO)  "  MAT0NEG(-7 , zl)  *  MAT0NEG(-13 , z2)  ; 
state_i  =  (state_i  +  31)  &  OxOOOOOOlfU; 
return  STATE [state_i] ; 

} 
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A.5  Multiply  With  Carry 
Complementary  Multiply  With  Carry  cmwc 


static  unsigned  int  Q[4096] ;  //  32  bit 

static  unsigned  long  long  c=362436;  //  64  bit 

unsigned  cmwc4096(  void  )  { 

//  static  unsigned  long  int  Q[4096] ; 

//  static  unsigned  long  int  c  =  362436; 
unsigned  long  long  int  t,  a=  1878211; 
static  unsigned  int  i  =  4095; 
unsigned  int  x,  r=  Oxfffffffe; 

i=(i+l)&  4095; 
t  =  a  *  Q[i]  +  c; 
c  =  (  t  »  32  ); 
x  =  t  +  c; 
if  (  x  <  c  )  { 
x++; 

C++; 

} 

return  (  Q[i]  =  r  -  x  ); 


Multiply  With  Carry  mwc 

static  unsigned  int  Q [2 56] ;  //  32  bit 

static  unsigned  long  long  carry=362436;  //  64  bit 

unsigned  MWC256(void)  { 

unsigned  long  long  t,a=809430660LL; 
static  unsigned  char  i=255; 

t=a*Q [++i] +carry ; 
carry=(t»32) ; 
return(Q[i]=t) ; 
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Multiply  With  Carry  mwcx 

static  unsigned  int  Q[256] ;  //  32  bit 

static  unsigned  long  long  c=362436;  //  64  bit 

unsigned  mwc256(  void  )  { 

unsigned  long  long  int  t,  a  =  154031582611 ; 

unsigned  int  x; 

static  unsigned  char  i  =  255; 

t  =  a  *  Q[++i]  +  c; 
c  =  (  t  »  32  ); 
x  =  t  +  c; 
if  (  x  <  c  )  { 
x++; 

C++; 

} 

return  (  Q [i]  =  x  ); 
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A.6  Exclusive-Or  (XOR)  Shift 
XOR  Shift  xorl28 

unsigned  int  xo8128(  void  )  {  //  return  an  integer  on  [0 , 2 ^ 32- 1] 

t=(x/'(x«ll)); 
x  =  y; 
y  -  z; 
z  =  w; 

w  =  (  w  (  w  »  19  )  )  *  (  t  '  (t»8)); 
return  w  ; 


A.7  Combined 

Keep  It  Simple  Stupid  kiss 

static  unsigned  int  x  =  123456789; 
static  unsigned  int  y  =  362436000; 
static  unsigned  int  z  =  521288629; 
static  unsigned  int  c  =  7654321; 

unsigned  int  kiss(  void  )  { 

unsigned  long  long  int  t,  a  =  69876906911; 

x  =  69069  *  x  +  12345; 

y  ~=  (  y  «  13  ) ; 

y  =  C  y  »  17  ) ; 

y  "=  (  y  «  5  ); 

t  =  a  *  z  +  c; 

c  =  (  t  »  32  ); 

return  x+y+  (z=t); 

} 
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Intentionally  Lett  Blank. 
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Appendix  B.  F2-Linear  RNG  Polynomial  Computation  Code 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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Intentionally  Lett  Blank. 
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B.l  Makefile 

##  makefile 

CPP=g++ 

CF=-02 

all:  poly_demo.exe  poly_mt32.exe  poly_tll3.exe  poly_mt64.exe  poly_t258.exe 

poly_%.exe:  poly_gen.cpp  rng_poly_%.h 

$(CPP)  $(CF)  -DRNG_POLY_H=\"rng_poly_$* . h\"  $<  -o  $@ 
strip  $@ 

clean: 

rm  -f  *.exe  * . exe . stackdump  *.bak  *“ 
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B.2  Driver  poly.gen .  cpp 


//  poly_gen. cpp 

#include  <iostream> 

#include  <i omani p> 

#include  <sstream> 

#include  <bitset> 

using  namespace  std; 

#include  RNG_POLY_H 

bool  cpgen  =  false;  //  generate  CP  flag 

int  bitnumber  =0;  //  bit  position  for  CP  calculation 

uz  seed  =  0xl2340f©0;  //  rng  seeed 

static  const  int  jbmax  =  20000;  //  binary  jump  representation  max  length 
int  jumplog  =20;  //  log_2  jump 

string  jumpstr  =  ;  //  jump  bit-string 

unsigned  long  long  jumpint  =  0;  //  jump  integer 

const  int  ubits  =  8*sizeof (uz) ;  //  bits  per  uz  (integer) 
const  int  wid  =  ubits/4;  //  uz  format  width 

string  USE  = 

"command  line  arguments :\n" 

"  CP  (characteristic  polynomial)  computation :\n" 

"  -c  :  compute  CP\n" 

"  -b  <n>  :  use  bit  position  n  to  generate  CP,  0  <=  n  <=  31  (or  63) \n" 

"  -s  <seed>  :  rng  seed,  hex  (0x...)»  for  CP  computation\n" 

"  JP  (jump  polynomial)  computation :\n" 

"  -js  <s>  :  jump  =  binary  string  s\n" 

"  -jl  <n>  :  jump  =  2/'n\n" 

"  -jn  <n>  :  jump  =  n  ,  decimal\n" 

"  -jx  <n>  :  jump  =  n  ,  hex\n"  ; 

//  Berlekamp-Massey  Algorithm 

void  BMA(  bitset<rng : :nS+l>&  c,  int&  L,  const  bitset<rng : :nS>&  s  ); 

//  jump  polynomial 

void  jpgen(  bitsetcrng: :CPx>&  jp,  const  bitset<jbmax>&  jumpbit,  const  bitset<rng : :CPx>&  r,  const  int&  n  ); 

//  set  jump  bits 
bitset<jbmax>  setjb(  ); 

//  format  bitset  to  hex  integers 

templatecint  N>  string  hexs(  const  bitset<N>&  x  ); 

int  main  (  int  argc  ,  char*  argv[]  )  { 

//  process  command  line  arguments 
if  (  argc  >  1  )  { 
stringstream  ss(""); 

for  (  char**  p  =  argv+1;  *p;  )  ss  «  "  "  «  *p++; 
string  s; 

while  (  !  ss.eofO  )  { 
ss  »  s; 

if  ("-c"  ==  s)  cpgen  =  true; 

else  if  ("-b"  ==  s)  ss  »  dec  »  bitnumber; 

else  if  ("-s"  ==  s)  ss  »  hex  »  seed; 

else  if  ("-jl"  ==  s)  {  //  integer  log_2  jump 
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ss  »  dec  »  jumplog,  jumpint  =  0,  jumpstr  = 

} 

else  if  C'-jn"  ==  s)  {  //  decimal  integer  jump 

jumplog  =  -1,  ss  »  dec  »  jumpint,  jumpstr  =  ; 

} 

else  if  C'-jx"  ==  s)  {  //  hex  integer  jump 

jumplog  =  -1,  ss  »  hex  »  jumpint,  jumpstr  = 

} 

else  if  C'-js"  ==  s)  {  //  bit  string  MSB...LSB  jump 

jumplog  =  -1,  jumpint  =0,  ss  »  jumpstr; 

} 

else  {  cout  «  "?  "  «  s  «  endl  «  USE;  exit(l);  } 

} 

} 

cout  «  "parameters:"  «  endl; 

if  (  cpgen  )  { 

cout  «  "seed  =  Ox"  «  setw(wid)  «  setfill(’©’)  «  hex  «  seed  «  endl 
«  "bitnumber  =  "  «  dec  «  bitnumber  «  endl  «  endl; 


else  { 

if  (  jumpstr  !=  ""  )  cout  «  "jumpstr  =  "  «  jumpstr  «  endl; 

if  (  jumplog  >  -1  )  cout  «  "jumplog  =  "  «  jumplog  «  endl; 

if  (  jumpint  )  cout  «  "jumpint  =  "  «  dec  «  jumpint  «  "d  =  0x"  «  setw(16) 

«  setfill(’Q’)  «  hex  «  jumpint  «  endl  «  endl; 


} 


//  compute  characteristic  polynomial 
if  (  cpgen  )  { 

cout  «  "compute  CP:"  «  endl; 

rng  x(  seed  ) ;  //  random  number  generator 

//  collect  bitstreams  at  position  bitnumber 
bitsetcrng : :nS>  b[rng::nG];  //  bitstream 
for  (  int  i=0;  i<rng::nS;  i++  ,  x.genC)  ) 
for  (  int  j=0;  j<rng::nG;  j++  ) 
b[j][i]  =  x.getstateO  [j]  »  bitnumber  &  1; 

//  compute  CP 
bitsetcrng : :nS+l>  CP; 
int  degCP; 

for  (  int  i=0;  i<rng::nG;  i++)  { 

BMA(  CP,  degCP,  b[i]  );  //  Berlekamp-Massey  Algorithm 

cout  «  "deg  =  "  «  degCP  «  endl  «  hexscrng : :nS+l>(CP)  «  endl; 

} 

return  0; 

} 


//  compute  jump  polynomial 

cout  «  "compute  JP"  «  endl  «  endl; 

bitset<jbmax>  jumpbit=setjb() ;  //  binary  jump  representation 
bitsetcrng: :CPx>  jp;  //  jump  polynomial 

for  (  int  i=©  ;  i  c  rng::nG;  i++)  { 

jpgenC  jp,  jumpbit,  rng::cp(i),  rng::cpd(i)  ); 
cout  cc  hexscrng : :CPx>(jp)  cc  endl  cc  endl; 

} 

}  //  end  of  mainQ 


void  BMA(  bitsetcrng : :nS+l>&  c,  int  &L,  const  bitsetcrng : :nS>&  s  )  { 
//  Berlekamp-Massey  Algorithm 

//  compute  the  minimal  polynomial  of  a  linearly  recurrent  sequence 

// 

//  LFSR  Synthesis  Algorithm  (Berlekamp  Iterative  Algorithm)  from 

// 
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//  J.  L.  Massey 

//  Shift-Register  Synthesis  and  BCH  Decoding 

// 

//  IEEE  Transactions  on  Information  Theory 

//  Volume  IT- 15  ,  Number  1  ,  January  1969  ,  pp  122-127 

// 

//  computes:  c  =  CP  and  L  =  linear  complexity  =  degree(c) 
int  d,  n  =  0,  x  =  1; 
bitsetcrng : :nS+l>  b(0),  t(0); 

L  =  0; 
c.resetO ; 
c[0]  =  b [0]  =  1; 
while  (  n  <  rng: :nS  )  { 
d  =  s  [n] ; 

for  (  int  i  =  1;  i  <=  L;  i++  ) 
d  *=  c[i]  &  s[n-i] ; 
if  C  !  d  ) 
x++; 

else  if  (  (  L«1  )  >  n  ) 
c  *=  b  «  x++; 
else  { 
t  =  c; 

c  *=  b  «  x  ; 

L  =  n  +  1  -  L; 
b  =  t; 
x  =  1; 

} 

n++; 

} 

//  linear  recurrence  reverses  massey’s  coefficients 
for  (  int  i=©  ;  i<=L/2 ;  i++  ) 
x=c[i] ,  c[i]=c[L-i],  c[L-i]=x; 


void  jpgen(  bitset<rng: :CPx>&  jp,  const  bitset<jbmax>&  jumpbit,  const  bitset<rng : :CPx>&  c,  const  int&  n  )  { 
//  computes:  jp(z)  =  jump  polynomialC  jump=jumpbit  ,  CP=c  ,  deg(CP)=  n  ) 


t; 


jp.resetO ; 
jp[l]=l; 

bitset<rng: :CPx> 
int  x; 

int  i  =  jbmax-1; 
while  (  ! (jumpbit) [i--]  ) 
for  (  ;  i>=0;  i —  )  { 
t .reset() ; 

for  (  int  j=n- 1 ;  j>=0; 
x  =  t [n-1] ; 
t  «=  1; 

if  (  x  )  t  *=  c; 
if  (  jp[j]  )  t  ~=  jp; 

} 

JP  =  t; 

if  (  (jumpbit) [i]  )  { 
x  =  jp[n-l] ; 
jp  «=  1; 
if  (  x  )  jp 

} 

} 

for  (  int  i=n;  i 
jp.reset(i) ; 


J” 


c; 


//  jp(z)  =  0 
//  jp(z)  =  z 
//  temporary  jp*jp 
//  carry  flag 
//  start  with  MSB 
//  skip  MSB 

//  apply  square-multiply  operations  to  jp 
//  square  jp 
)  { 


//  jp(z)  =  jp(z)  *  jp(z) 
//  multiply  jp(z)  by  z 

//  jp(z)  =  jp(z)  *  z 


<  rng::CPx;  i++)  //  zero  extra  high  bits 
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bitset<jbmax>  setjb(  )  { 

//  set  binary  expansion  of  jump  for  square-multiply  sequence 
if  (  jumpint  >  0  )  //  jump  =  unsigned  long  long  jump 

return  bitset<jbmax>  (jumpint); 
if  (  jumpstr  !=  ""  )  //  jump  =  bit  string  jumpstr 

return  bitset<jbmax>  (jumpstr); 
bitset<jbmax>  j;  //  jump  =  Z^jumplog 

j . set(jumplog) ; 
return  j ; 

} 

template<int  N>  string  hexs(  const  bitset<N>&  x  )  { 
stringstream  ss(""); 
uz  z,  m; 

int  j,  jj,  k,  n=64/wid; 

for  (  j=jj=0;  j  <  rng::nW;  j++  )  { 

for  (  z=k=0,  m=l;  k<ubits;  k++,  jj++,  m«=l) 
if  (  jj  <  x.sizeO  &&  x[jj]  ) 
z  |=  m; 

ss  «  "Ox"  «  hex  «  setw(wid)  «  setfill  (’O’)  «  z; 
if  (  j  %  n  ==  n-1  )  ss  «  endl; 
else  if  (  j  <  rng::nW-l  )  ss  «  "  "; 

//  else  ss  «  "  "; 

1 

return  ss.strQ; 
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B.3  Demonstration  demo  RNG  class  rng_poly_demo .  h 

//  -"-C++-"- 

//  rng_poly_derao .h 

typedef  unsigned  int  uz; 

class  rng  { 
public: 

rng(  uz  seed®  )  ; 
uz  gen(  )  ; 
uz*  getstateC  ); 
static  const  int  nG  =  1; 
static  const  int  nW  =  1; 
static  const  int  nS  =  100; 
static  const  int  CPx  =  31; 
static  int  cpd(  int  )  ; 
static  bitset<CPx>  cp(  int 
private : 

struct  state  { 
uz  z  [1] ; 

}; 

state  s; 

void  seed(  uz  s®  ) ; 
static  const  int  CPd; 
static  const  uz  CP; 

}; 

rng:: rng (  uz  seed®  )  { 
seed(  seed®  ) ; 
gen() ; 

} 

void  rng::seed(  uz  s®  )  { 
s.z[Q]  =  sQ; 

} 

//  Tausworthe  t®88.®  (the  first  generator),  period  =  2 "3 1-1 
uz  rng : : gen(  )  { 

*s.z  =  ((*s.z  &  Qxfffffffe)  «  12  )  "  (  (  (*s.z  «  13)  "  *s.z  )  »  19) 
return  *s.z; 

} 

uz*  rng :: getstateC  )  {  //  return  state 

return  s.z; 

} 

bitsetcrng : :CPx>  rng::cp(  int  )  {  //  return  CP  bitset 
return  bitsetcrng: :CPx>  (  CP  ) ; 

} 

int  rng::cpd(  int  )  {  //  return  CP  degree 

return  CPd; 

} 

//  characteristic  polynomial  degree 
const  int  rng:: CPd  =31  ; 

//  characteristic  polynomial  coefficients 
const  uz  rng::CP  =  ®x82Q820®lu; 


//  constructor 
//  generator 
//  return  state 

//  number  of  generators  (  CPs  ) 

//  words  in  a  generator  (  CP  ) 

//  number  of  states  for  CP  computation 
//  CP  max  degree 
//  return  CP  degree 
) ;  //  return  CP 

//  rng  state  structure 


//  rng  state 
//  set  seed 
//  CP  degree 
//  CP 


B.4  Tausworthe  tll3  RNG  class  rng_poly_tll3  .h 


//  -*-C++-*- 

//  rng_poly_tll3 .h 

typedef  unsigned  int  uz; 

class  rng  { 
public: 

rng(  uz  seed®  );  //  constructor 

uz  gen(  ) ;  //  generator 

uz*  getstateC  );  //  return  state 

static  const  int  nG  =  4;  //  number  of  generators  (  CPs  ) 

static  const  int  nW  =  1;  //  words  in  a  generator  (  CP  ) 

static  const  int  nS  =  100 ;  //  number  of  states  for  CP  computation 

static  const  int  CPx  =31;  //CP  max  degree 

static  int  cpd(  const  int&  n  ) ;  //  return  CP  degree 
static  bitset<CPx>  rng::cp(  const  int&  n  ) ;  //  return  CP 
private : 

struct  state  {  //  rng  state  structure 

uz  z [4] ; 

}; 

state  s;  //  rng  state 

void  seed(  uz  s®  ) ;  //  set  seed 

static  const  uz  CP[nG];  //  CP 

static  const  int  CPd[nG] ;  //  CP  degree 

}; 

rng:: rng (  uz  seed®  )  { 
seed(  seed®  ) ; 
gen() ; 


void  rng::seed(  uz  s®  )  { 
s.z[Q]  =  s®  *  69069; 
if  (  s.z[0]  <  2  )  s.z[Q]  +=  2U; 
s.z[l]  =  s.z[Q]  *  69069; 
if  (  s.z[l]  <  8  )  s.z[l]  +=  8U; 
s.z[2]  =  s.z[l]  *  69069; 
if  (  s.z[2]  <  16  )  s.z[2]  +=  16U; 
s.z[3]  =  s.z[2]  *  69069; 
if  (  s.z[3]  <  128  )  s.z[3]  +=  128U; 


uz  rng: :gen(  )  { 

s.z[Q]  =  ((s.z[Q]  &  Qxfffffffe)  «  18)  "  (((s.z[Q]  «  6)  *  s.z[Q])  »  13) 

s.z[l]  =  C(s.z[l]  &  ®xfffffff8)  «  2)  ~  (C(s.z[l]  «  2)  *  s.z[l])  »  27) 

s.z[2]  =  C(s.z[2]  &  Oxfffffff®)  «  7)  ~  C((s.z[2]  «  13)  "  s.z[2])  »  21) 

s.z[3]  =  C(s.z[3]  &  ®xffffff80)  «  13)  ~  C((s.z[3]  «  3)  ~  s.z[3])  »  12) 

return  (  s.z[®]  "  s.z[l]  *  s.z[2]/'  s.z[3]  )  ; 

} 

uz*  rng :: getstateC  )  {  //  return  state 

return  s.z; 

} 

bitsetcrng : :CPx>  rng::cp(  const  int&  n  )  {  //  return  CP  bitset 
return  bitsetcrng: :CPx>  (  CP[n]  ); 

} 

int  rng::cpd(  const  int&  n  )  {  //  return  CP  degree 

return  CPd[n] ; 
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} 

//  characteristic  polynomial  degree 

const  int  rng::CPd[nG]  =  {  31  ,  29  ,  28  ,  25  }  ; 

//  characteristic  polynomial  coefficients 
const  uz  rng : : CP [nG]  =  { 

©x80400855u  ,  0x2Q00Q005u  ,  0xllll3111u  ,  0x02041879u 

}  ; 
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B.5  Tausworthe  t2 58  RNG  class  rng_poly_t2 58 . h 


//  -*-C++-*- 

//  rng_poly_t258 .h 


typedef  unsigned  long  long  int  uz; 


class  rng  { 
public: 

rng(  uz  seed®  ); 
uz  gen(  ); 
uz*  getstateC  ); 
static  const  int  nG 
static  const  int  nW 
static  const  int  nS 


//  constructor 
//  generator 
//  return  state 

5;  //  number  of  generators  (  CPs  ) 

1;  //  words  in  a  generator  (  CP  ) 

200;  //  number  of  states  for  CP  computation 

//  CP  max  degree 


static  const  int  CPx  =  63; 
static  int  cpd(  const  int&  n  ) ;  //  return  CP  degree 
static  bitset<CPx>  cp(  const  int&  n  ) ;  //  return  CP 
private : 

struct  state  {  //  rng  state  structure 

uz  z[5] ; 


}; 

state  s; 

void  seed(  uz  s®  ) ; 
static  const  int  CPd[nG] ; 
static  const  uz  CP[nG]; 


//  rng  state 
//  set  seed 
//  CP  degree 
//  CP 


rng:: rng (  uz  seed®  )  { 
seed(  seed®  ) ; 
gen() ; 

} 


void  rng::seed(  uz  s®  )  { 
s.z[Q]  =  s®  *  69069; 
if  (  s.z[Q]  <  2  )  s.z[Q]  +=  2u; 
s.z[l]  =  s.z[Q]  *  69069; 
if  (  s.z[l]  <  8  )  s.z[l]  +=  8u; 
s.z[2]  =  s.z[l]  *  69069; 
if  (  s.z[2]  <  16  )  s.z[2]  +=  16u; 
s.z[3]  =  s.z[2]  *  69069; 
if  (  s.z[3]  <  128  )  s.z[3]  +=  128u; 
s.z[4]  =  s.z[3]  *  69069; 
if  (  s.z[4]  <  2048  )  s.z[4]  +=  2®48u; 

} 


uz  rng : : gen(  )  { 


s.z[®]  =  C(s.z[0] 

&  Oxfffffffffffffffeull) 

« 

is)  ■ 

(C(s.z[0] 

« 

i) 

s.z[0]) 

» 

53) 

s . z  [  1]  =  C(s.z[l] 

&  OxfffffffffffffeOOull) 

« 

5)  ' 

C((s.z[l] 

« 

24) 

■  s.z[l]) 

» 

50) 

s.z[2]  =  C(s.z[2] 

&  QxfffffffffffffQQQull) 

« 

29)  " 

(((s.z [2] 

« 

3) 

-  s . z [2] ) 

» 

23) 

s.z  [3]  =  C(s.z[3] 

&  OxfffffffffffeOOOOull) 

« 

23)  - 

(  ( (s . z [3] 

« 

5) 

"  s . z [3] ) 

» 

24) 

s.z[4]  =  C(s.z[4] 

&  ®xffffffffff8®0®®®ull) 

« 

8)  ' 

(((s.z [4] 

« 

3) 

s.z [4]) 

» 

33) 

return  (  s.z[®] 

s.z[l]  ~  s.z[2]  "  s.z[3] 

:.z[4] 

)  ; 

} 

uz*  rng :: getstateC  )  {  //  return  state 

return  s.z; 

} 

bitsetcrng : :CPx>  rng::cp(  const  int&  n  )  {  //  return  CP  bitset 
bitsetcrng : :CPx>  c  =  Q; 

(  c  |=  CP[n]  »  32  )  «=  32; 
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C  1=  CP [n] ; 
return  c; 

} 

int  rng: :cpd(  const  int&  n  )  {  //  return  CP  degree 
return  CPd[n] ; 

} 

//  characteristic  polynomial  degree 

const  int  rng::CPd[nG]  =  {  63  ,  55  ,  52  ,  47  ,  41  }  ; 

//  characteristic  polynomial  coefficients 
const  uz  rng : : CP [nG]  =  { 

0X8000004O00002003U11  ,  0X0080100001000801U11  ,  ©x®01®O©08®805414dull  , 
0x®00®8024092248blull  ,  0X000002OO00O00O09U11 

}; 
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B.6  Mersenne  Twister  mt32  RNG  class  rng_polyjtit32  .h 


//  -*-C++-*- 
//  rng_poly_mt32 .h 

typedef  unsigned  int  uz; 

class  rng  { 
public: 

rng(  uz  seed®  );  //  constructor 

uz  gen(  ) ;  //  generator 

uz*  getstateC  );  //  return  state 

static  const  int  nG  =  1;  //  number  of  generators  (  CPs  ) 

static  const  int  nW  =  624;  //  words  in  a  generator  (  CP  ) 

static  const  int  nS  =  5000® ;  //  number  of  states  for  CP  computation 

static  const  int  CPx  =  19937;  //  CP  max  degree 

static  int  cpd(  int  ) ;  //  return  CP  degree 

static  bitset<CPx>  cp(  int  );  //  return  CP 
private : 

struct  state  {  //  rng  state  structure 

int  n; 
uz  y [1] ; 
uz  z [624] ; 

}; 

state  s;  //  rng  state 

void  seed(  uz  seed®  ) ;  //  set  seed 

static  const  uz  CP[nW];  //  CP 

static  const  int  CPd;  //  CP  degree 

static  const  uz  N  =  624;  //  rng  parameter 

static  const  uz  M  =  397 ;  //  rng  parameter 

static  const  uz  A  =  Qx99Q8bQdf;  //  constant  vector  a 

static  const  uz  UM  =  0x8000000®;  //  most  significant  w-r  bits 

static  const  uz  LM  =  0x7fffffff;  //  least  significant  r  bits 

}; 

rng:: rng (  uz  seed®  )  { 
seed(  seed®  ) ; 
gen() ; 

} 

void  rng::seed(  uz  seed®  )  { 
s.z[Q]=  seed®; 

for  (  s.n=l  ;  s.n  <  N  ;  s.n++  ) 

s.z[s.n]  =  1812433253  *  (s.z[s.n-l]  "  (s.z[s.n-l]  »  3®))  +  s.n; 

} 

uz  rng: :gen(  )  { 
uz  y; 

static  uz  magOl [2]={®x® ,  A}; 

if  (s.n  >=  N)  {  //  generate  N  words  at  one  time 

int  kk; 

for  (  kk=Q  ;  kk<N-M  ;  kk++  )  { 

*s.y  =  (  s.z[kk]  &  UM  )  |  (  s.z[kk+l]  &  LM  ) ; 
s.z[kk]  =  s.z[kk+M]  "  (*s.y  »  1)  ~  mag®l[*s.y  &  ®xl]; 

} 

for  (  ;  kk<N-l  ;  kk++  )  { 

*s.y  =  (  s.z[kk]  &  UM  )  |  (  s.z[kk+l]  &  LM  ) ; 

s.z[kk]  =  s.z[kk+(M-N)]  ~  (*s.y  »  1)  *  mag®l[*s.y  &  ®xl]; 

} 

*s.y  =  (  s.z[N-l]  &  UM  )  |  (  s.z[Q]  &  LM  ); 
s.z[N-l]  =  s.z[M-l]  *  (*s.y  »  1)  ~  mag®l[*s.y  &  ®xl]; 
s.n  =  Q; 
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} 

y  =  *s.y  =  s.z[s.n++]; 

//  Tempering 
y  “=  (y  »  11); 
y  ~  =  (y  «  7)  &  0x9d2c5680; 

y  ~=  (y  «  15)  &  0xefc600©0; 
y  ~=  (y  »  18) ; 
return  y; 


uz*  rng : :getstate(  )  {  //  return  state 

return  s.y; 

} 

bitsetcrng : :CPx>  rng::cp(  int  )  {  //  return  CP  bitset 
bitsetcrng : :CPx>  c  =  0; 
for  (  int  j=nW-l;  j;  j — ) 

C  c  |=  CP[j]  )  «=  32; 
c  |=  CP[0] ; 
return  c; 

} 

int  rng::cpd(  int  )  {  //  return  CP  degree 

return  CPd; 

} 

//  characteristic  polynomial  degree 

const  int  rng:: CPd  =  19937; 
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//  characteristic  polynomial  coefficients 
const  uz  rng : : CP [nW]  =  { 
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B.7  Mersenne  Twister  mt64  RNG  class  rng.poly jnt64 .  h 


//  -*-C++-*- 
//  rng_poly_mt64 .h 

typedef  unsigned  long  long  int  uz; 

class  rng  { 
public: 

rng(  uz  seed®  );  //  constructor 

uz  gen  (  ) ;  //  generator 

uz*  getstateC  );  //  return  state 

static  const  int  nG  =  1;  //  number  of  generators  (  CPs  ) 

static  const  int  nW  =  312;  //  words  in  a  generator  (  CP  ) 

static  const  int  nS  =  50000;  //  number  of  states  for  CP  computation 

static  const  int  CPx  =  19937;  //  CP  max  degree 

static  int  cpd(  int  ) ;  //  return  CP  degree 

static  bitset<CPx>  cp(  int  );  //  return  CP 
private : 

struct  state  {  //  rng  state  structure 

int  n; 
uz  y [1] ; 
uz  z  [312] ; 

}; 

state  s;  //  rng  state 

void  seed(  uz  seed©  ) ;  //  set  seed 

static  const  uz  CP[nW];  //  CP 

static  const  int  CPd;  //  CP  degree 

static  const  uz  NN  =  312;  //  rng  parameter 

static  const  uz  MM  =  156;  //  rng  parameter 

static  const  uz  MATRIX_A  =  0xB5026F5AA96619E9ull ;  //  constant  vector  a 

static  const  uz  UM  =  0xFFFFFFFF80000000ull ;  //  most  significant  w-r  bits 

static  const  uz  LM  =  0x000000007FFFFFFFull ;  //  least  significant  r  bits 

}; 

rng:: rng (  uz  seed®  )  { 
seed(  seed®  ) ; 
gen() ; 

} 

void  rng::seed(  uz  seed®  )  { 
s.z[0]  =  seedO; 

for  (  s.n  =  1  ;  s.n<NN  ;  s.n++) 

s.z[s.n]  =  6364136223846793005ull  *  (s.z[s.n-l]  *  (s.z[s.n-l]  »  62))  +  s.n; 

} 

uz  rng : : gen  (  ) 

{ 

uz  y; 

static  uz  magOl [2]={®x0ull ,  MATRIX_A} ; 

if  (s.n  >=  NN)  {  //  generate  N  words  at  one  time 

int  kk; 

for  (  kk  =  0  ;  kk  <  NN  -  MM  ;  kk++  )  { 

*s.y  =  (  s.z[kk]  &  UM  )  |  (  s.z[kk+l]  &  LM  ) ; 

s.z[kk]  =  s.z[kk+MM]  "  (*s.y  »  1)  "  mag01[(int)  (*s.y  &  Oxlull)]; 

} 

for  (  ;  kk  <  NN  -  1  ;  kk++  )  { 

*s.y  =  (  s.z[kk]  &  UM  )  |  (  s.z[kk+l]  &  LM  ) ; 

s.z[kk]  =  s. z[kk+(MM-NN)]  "  (*s.y  »  1)  ~  mag01[(int)  (*s.y  &  Oxlull)]; 

} 

*s.y  =  (  s.z[NN-l]  &  UM  )  |  (  s.z[0]  &  LM  ) ; 

s.z[NN-l]  =  s.z[MM-l]  "  (*s.y  »  1)  "  mag01[  (int)  (*s.y  &  Oxlull)]; 
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s.n  =  0; 

} 

y  =  *s.y  =  s.z[s.n++]; 

//  Tempering 

y  ~=  (y  »  29)  &  Ox5555555555555555ull 
y  ~  =  (y  «  17)  &  Qx71D67FFFEDA6QQQQull 
y  ~=  (y  «  37)  &  0xFFF7EEE000000000ull 
y  ~=  (y  »  43)  ; 
return  y; 


uz*  rng : :getstate(  )  {  //  return  state 

return  s.y; 

} 


bitsetcrng : :CPx>  rng::cp(  int  )  {  //  return  CP  bitset 
bitsetcrng: :CPx>  c  =  0; 
for  (  int  j=nW-l;  j;  j — )  { 

(  c  |=  CP[j]  »  32  )  «=  32; 

(  c  |=  CP [ j ]  )  «=  32; 

} 

(  c  |=  CP[Q]  »  32  )  «=  32; 
c  |=  CP[0] ; 
return  c; 


int  rng::cpd(  int  )  {  //  return  CP  degree 

return  CPd; 

} 


//  characteristic  polynomial  degree 
const  int  rng:: CPd  =  19937; 
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//  characteristic  polynomial  coefficients 
const  uz  rng : : CP [nW]  =  { 


QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 


QxQQQQQQQQQQQQQQQQull 

QxQQlQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxSQQQQQQQSSSQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQlull 

QxQQ8QQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

Qx4QQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQlQ58QQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQlQQ886QQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQlQQ5QQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQlQQQ8ull 

QxQQQ4QQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQlQull 

QxQQ4QQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQ4QQ4QQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQ4QQQ4Q4QQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQ4ull 

QxQQQ8Q224QQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQ4ull 

QxQ4QQ4QQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQ82624ull 

Qx84QQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQ4Q4QQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQ184QQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQ4Q4QQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQ5ull 

QxQ4a4QQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQ4Qull 

Qx4QQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQa4QQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQ4QQQQQQQull 


QxQQQQQQQQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 

QxQQQQQ4QQQQQQQQQQull 

QxQQQQQQQQQQQQQQQQull 


}  ; 
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Appendix  C.  Independent  F2-Linear  RNG  Implementation 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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Intentionally  Lett  Blank. 
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C.l  Makefile 

##  makefile 

CF  =  -02 
CPP  =  g++ 

all:  jump_demo.exe  jump_tll3.exe  jump_mt32.exe  jump_t258.exe  jump_mt64.exe 

j  ump_% . exe :  j  ump_gen . cpp  rng_ j ump_% . h 

$(CPP)  $(CF)  -DRNG_JUMP_H=\"rng_jump_$* .h\"  $<  -o  $@ 
strip  $@ 

clean: 

-rm  -f  *.o  *.exe 


53 


C.2  Driver  jump _gen.cpp 


//  -*-C++-*- 

//  jump_gen. cpp 
//  driver  for  rng  class 

//  Random  Number  Generator  with  statistically  independent  instances 

#include  <iostream> 

#include  <sstream> 

#include  <iomanip> 

using  namespace  std; 

#include  RNG_JUMP_H 

#define  USE()  cout  «  "args  are:  -s  seed  -n  number_of_rngs"  «  endl 

int  main(  int  argc  ,  char*  argv[]  )  { 

int  n_gen  =5;  //  number  of  independent  rngs 

//  default  seed 
unsigned  seed=Oxl2340f©0 ; 

if  (  argc  >  1  )  { 
stringstream  ss(""); 

for  (  char**  p=argv+l  ;  *p  ;  )  ss  «  "  "  «  *p++; 
string  tok; 

while  (  !  ss.eofC)  )  { 
ss  »  tok; 

if  C  "-n"  ==  tok  )  ss  »  dec  »  n_gen; 

else  if  (  "-s"  ==  tok  )  ss  »  hex  »  seed; 

else  {  cout  «  "?  "  «  tok  «  endl  «  USE();  exit(l);  } 

} 

} 

//  set  the  seed  before  constructing  rngs 
rng : : setseedC  seed  ); 

//  create  rngs 
rng  Z [n_gen] ; 

//  call  rngs  and  print 

for  (  int  i=©  ;  i  <  n_gen  ;  i++  ) 

cout  «  "call  Z["  «  dec  «  setw(2)  «  setfill  (’  ’)  «  i  «  "].gen() 
«  hex  «  setw(rng : :wid)  «  setfill  (’0’)  «  Z[i].gen()  «  endl; 

//  advance  rngs  to  the  same  state 
cout  «  "synchronize  the  rngs"  «  endl; 
for  (  int  i=©  ;  i  <  n_gen  ;  i++  ) 

for  (  int  j  =  0  ;  j  <  (n_gen-i-l)*(l«20)  ;  j++  ) 

Z[i] .gen() ; 

//  call  rngs  and  print 

for  (  int  i=©  ;  i  <  n_gen  ;  i++  ) 

cout  «  "call  Z["  «  dec  «  setw(2)  «  setfill  (’  ’)  «  i  «  "].gen() 
«  hex  «  setw(rng : :wid)  «  setfill  (’O’)  «  Z[i].gen()  «  endl; 

return  0; 

} 
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C.3  Demonstration  demo  RNG  class  rng_jump_demo . h 


//  -*-C++-*- 

//  rng_ j ump_derao . h 

//  F2-Linear  Feedback  Shift  Register  Random  Number  Generator 
//  with  statistically  independent  instances 
#include  <bitset> 

class  rng  { 

struct  state  {  //  state  information 

unsigned  int  z; 

}  I 

public: 

//  the  first  instance  is  initialized  from  the  seed 
//  subsequent  instances  are  jumped  ahead  giving  independent  rngs 
rng(  );  //  constructor  (called  with  no  arguments) 

rng(  const  unsigned&  seed  );  //  constructor  (called  with  seed  argument) 

unsigned  gen(  );  //  rng  returns  a  32-bit  unsigned  integer 

static  void  setseed(  const  unsigned&  seed  ) ;  //  call  once  for  all  rngs 
static  const  int  wid  =8;  //  integer  format  width 

private : 

state  s  ;  //  rng  state 

unsigned  id  ;  //  rng  sequential  id  per  instance 

void  init(  );  //  initialize  a  rng 

void  seedgen(  ) ;  //  seed  the  first  rng 

void  jump(  );  //  jump  to  new  state 

void  jpunpack(  const  unsigned&  Jump_P  ); 

static  unsigned  seed;  //  seed  for  first  generator 

static  unsigned  ref_counter;  //  number  of  instances 

static  state  ref_state;  //  zero  state  of  newest  instance 

static  unsigned  Jump_P[];  //  jump  polynomial  coefficients 

static  bitset<32>  JP;  //  jump  polynomial 

}; 

rng:: rng (  )  {  //  constructor  (called  with  no  arguments) 

init(  ); 

} 

rng:: rng (  const  unsigned&  seed  )  {  //  constructor  (called  with  seed  argument) 
setseed(  seed  ) ; 
init(  ); 

} 

void  rng : : setseed(  const  unsigned&  seed  )  { 
if  (  ref_counter  ==  0  ) 
rng : : seed  =  seed; 

} 

unsigned  rng::gen(  )  {  //  rng  returns  a  32-bit  unsigned  integer 

//  this  is  Tausworthe  t088.0  (the  first  generator),  period  =  2^31-1 
s.z  =  (  (s.z  &  Oxfffffffe)  «  12  )  ~  (  (  (s.z  «  13)  *  s.z  )  »  19  ); 
return  s.z; 

} 

void  rng : : init(  )  { 

//  the  first  rng  state  s  is  seeded  as  indicated 
//  for  subsequent  rngs  compute  state  s 
//  jumped  ahead  from  the  reference  state 

//  set  reference  state  for  next  jump  to  zero  state  of  current  rng 
if  (  ref_counter  ==  ©  )  { 
seedgen(  ) ; 

jpunpack(  Jump_P[0]  ); 
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} 

else 
jump(  ); 
ref_state  =  s  ; 
id  =  ref_counter++; 


void  rng : : seedgenC  )  { 

//  seed  the  system 

//  set  seed  for  first  generator  (algorithm  specific  to  rng) 
//  the  rng  must  be  in  state  zero,  not  the  seed  state 
s.z  =  seed; 
gen(  ); 


void  rng : : jpunpack(  const  unsigned&  Jump_P  )  { 

JP  |=  Jump_P; 

} 

void  rng::jump(  )  { 

//  set  rng  state  s  to  reference  state 

//  compute  temp_state  =  J[0]*s[0]  +  ...  +  J[31]*s[31] 

//  set  s  to  temp_state 

//  now  s  is  jumped  ahead  from  reference  state  by  J 
state  temp_state; 
s  =  ref_state; 
temp_state.z  =  0; 

for  (  int  i  =  0  ;  i  <  32  ;  i++  )  { 
if  (  JP[i]  ) 
temp_state.z  "=  s.z; 
gen(  ) ; 

} 

s  =  temp_state  ; 

} 

//  number  of  rngs  instantiated 
unsigned  rng: :ref_counter=Q; 

//  reference  state  for  jump  computation  is  the 
//  zero  state  of  the  previous  instance 
rng : : state  rng : : ref_state ; 

//  seed  for  first  rng  instance 
unsigned  rng : : seed=0xl2340f0Q ; 

//  Jump  polynomial 
bitset<32>  rng::JP  =  0; 

//  jump  polynomial  coefficients  for  jump  =  2 "20 
unsigned  rng : : Jump_P[l]  =  { 

0x4fbc5320 

}; 
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C.4  Tausworthe  tll3  RNG  class  rng_jump_tll3  .h 


//  -*-C++-*- 

//  rng_ j ump_t 1 1 3 . h 

//  Tausworthe  Random  Number  Generator 
//  with  statistically  independent  instances 
#include  <cstring> 

#include  <bitset> 

class  rng  { 

struct  state  {  //  state  information 

unsigned  int  z [4] ; 

}  I 

public: 

//  the  first  instance  is  initialized  from  the  seed 
//  subsequent  instances  are  jumped  ahead  giving  independent  rngs 
rng(  );  //  constructor  (called  with  no  arguments) 

rng(  const  unsigned&  seed  );  //  constructor  (called  with  seed  argument) 

unsigned  gen(  );  //  rng  returns  a  32-bit  unsigned  integer 

static  void  setseed(  const  unsigned&  seed  ) ;  //  call  once  for  all  rngs 
static  const  int  wid  =8;  //  integer  format  width 

private : 

state  s  ;  //  rng  state 

unsigned  id  ;  //  rng  sequential  id  per  instance 

void  init(  );  //  initialize  a  rng 

void  seedgen(  ) ;  //  seed  the  first  rng 

void  jump(  );  //  jump  to  new  state 

void  jpunpack(  const  unsigned*  const  Jump_P  ); 

static  unsigned  seed;  //  seed  for  first  generator 

static  unsigned  ref_counter;  //  number  of  instances 

static  state  ref_state;  //  zero  state  of  newest  instance 

static  unsigned  Jump_P[][4];  //  jump  polynomial  coefficients 

static  bitset<32>  JP[4];  //  jump  polynomial 

}; 


rng:: rng (  )  {  //  constructor  (called  with  no  arguments) 

init(  ); 

} 

rng:: rng (  const  unsigned&  seed  )  {  //  constructor  (called  with  seed  argument) 
setseed(  seed  ) ; 
init(  ); 

} 

void  rng : : setseed(  const  unsigned&  seed  )  { 
if  (  ref_counter  ==  0  ) 
rng : : seed  =  seed; 

} 

unsigned  rng::gen(  )  {  //  rng  returns  a  32-bit  unsigned  integer 

s.z[0]  =  ((s.z[0]  &  Oxfffffffe)  «  18)  "  (((s.z[0]  «  6)  "  s.z[0])  »  13); 

s . z [  1]  =  ((s.z[l]  &  0xfffffff8)  «  2)  ~  (((s.z[l]  «  2)  "  s . z [  1] )  »  27); 

s.z[2]  =  ((s.z[2]  &  OxfffffffO)  «  7)  "  (((s.z[2]  «  13)  "  s.z[2])  »  21); 

s.z[3]  =  ((s.z[3]  &  0xffffff80)  «  13)  "  (((s.z[3]  «  3)  *  s.z[3])  »  12); 

return  (  s.z[0]  "  s.z[l]  "  s.z[2]~  s.z[3]  )  ; 

} 

void  rng : : init(  )  { 

//  the  first  rng  state  s  is  seeded  as  indicated 
//  for  subsequent  rngs  compute  state  s 
//  jumped  ahead  from  the  reference  state 
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//  set  reference  state  for  next  jump  to  zero  state  of  current  rng 
if  (  ref_counter  ==  0  )  { 
seedgenC  ) ; 

jpunpack(  Jump_P[0]  ); 

} 

else 
jumpC  ); 
ref_state  =  s  ; 
id  =  ref_counter++; 


void  rng :: seedgenC  )  { 

//  seed  the  system 

//  set  seed  for  first  generator  (algorithm  specific  to  rng) 
//  the  rng  must  be  in  state  zero,  not  the  seed  state 


s.z[0]  =  seed 
if  (  s.z[0]  < 
s.z[l]  =  s.z[0] 
if  (  s.z[l]  < 
s.z[2]  =  s.z[l] 
if  (  s.z[2]  < 


*  69069; 

2  )  s.z[Q]  +=  2U; 

*  69069; 

8  )  s . z [ 1]  +=  8U; 

*  69069; 

16  )  s.z[2]  +=  16U; 


s.z[3]  =  s.z[2]  *  69069; 
if  (  s.z[3]  <  128  )  s . z[3] 
gen() ; 


h=  128U; 


void  rng : : jpunpackC  const  unsigned*  const  Jump_P  )  { 
for  (  int  j=0;  j  <4 ;  j++  ) 

JP[j]  |=Jump_P[j]; 

} 


void  rng:: jumpC  )  { 

//  set  rng  state  s  to  reference  state 

//  compute  temp_state  =  J[0]*s[0]  +  ...  +  J[31]*s[31] 

//  set  s  to  temp_state 

//  now  s  is  jumped  ahead  from  reference  state  by  J 
int  i  ,  j ; 
state  temp_state; 
s  =  ref_state; 

memsetC  temp_state.z  ,  0  ,  sizeof (temp_state . z)  ); 
for  (  i  =  0  ;  i  <  32  ;  i++  )  { 
for  (j=0  ;  j<4  ;  j++) 
if  (  JP[j][i]  ) 

temp_state.z[j]  ~=  s . z [ j ] ; 
gen() ; 

} 

s  =  temp_state  ; 


//  number  of  rngs  instantiated 
unsigned  rng: :ref_counter=0; 

//  reference  state  for  jump  computation  is  the 
//  zero  state  of  the  previous  instance 
rng:: state  rng: :ref_state; 

//  seed  for  first  rng  instance 
unsigned  rng : : seed=0xl2340f00 ; 

//  Jump  polynomial 

bitset<32>  rng::JP[4]  =  {  0,  0,  0,  0  }; 
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//  Tausll3  jump  coefficients 
unsigned  int  rng : : Jump_P [1] [4]  =  { 

{ 

//  2"2S 

0x®c382e31  ,  ®xlb®4®425  ,  8x®b49a5®9  ,  0xS173f6b® 

//  2"8S 

//  8x487cf69c  ,  ®x®8be631®  ,  8x®4bfe2bb  ,  ®x8®Q824f9 

} 

}; 
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C.5  Tausworthe  t2 58  RNG  class  rng_jump_t2 58 . h 


//  -*-C++-*- 

//  rng_jump_t258 .h 

//  Tausworthe  Random  Number  Generator 
//  with  statistically  independent  instances 
#include  <cstring> 

#include  <bitset> 

class  rng  { 

struct  state  {  //  state  information 

unsigned  long  long  int  z [5] ; 

}  I 

public: 

//  the  first  instance  is  initialized  from  the  seed 
//  subsequent  instances  are  jumped  ahead  giving  independent  rngs 
rng(  );  //  constructor  (called  with  no  arguments) 

rng(  const  unsigned  long  long  int&  seed  );  //  constructor  (called  with  seed  argument) 

unsigned  long  long  int  gen(  ) ;  //  rng  returns  a  64-bit  unsigned  integer 
static  void  setseed(  const  unsigned  long  long  int&  seed  ) ;  //  call  once  for  all  rngs 
static  const  int  wid  =16;  //  integer  format  width 

private : 

state  s  ;  //  rng  state 

unsigned  id  ;  //  rng  sequential  id  per  instance 

void  init(  );  //  initialize  a  rng 

void  seedgen(  ) ;  //  seed  the  first  rng 

void  jump(  );  //  jump  to  new  state 

void  jpunpack(  const  unsigned  long  long*  const  Jump_P  ); 

static  unsigned  long  long  int  seed;  //  seed  for  first  generator 

static  unsigned  ref_counter;  //  number  of  instances 

static  state  ref_state;  //  zero  state  of  newest  instance 

static  unsigned  long  long  Jump_P[] [5] ;  //  jump  polynomial  coefficients 

static  bitset<64>  JP [5] ;  //  jump  polynomial 

}; 

rng::rng(  )  {  //  constructor  (called  with  no  arguments) 

init(  ); 

} 

rng:: rng (  const  unsigned  long  long&  seed  )  {  //  constructor  (called  with  seed  argument) 
setseed(  seed  ) ; 
init(  ); 


void  rng : : setseed(  const  unsigned  long  long&  seed  )  { 
if  (  ref_counter  ==  0  ) 
rng : : seed  =  seed; 

} 

unsigned  long  long  rng::gen(  )  {  //  rng  returns  a  64-bit  unsigned  integer 

s.z[Q]  =  ((s.z[0]  &  Oxfffffffffffffffeull)  «  10)  ~  (((s.z[Q]  «  1)  ~  s.z[Q])  »  53) 

s . z [1]  =  ((s.z[l]  &  OxfffffffffffffeOOull)  «  5)  ~  (((s.z[l]  «  24)  *  s.z[l])  »  50) 

s.z[2]  =  ((s.z[2]  &  OxfffffffffffffOOOull)  «  29)  ~  (((s.z[2]  «  3)  *  s.z[2])  »  23) 

s.z[3]  =  ((s.z[3]  &  OxfffffffffffeOOOOull)  «  23)  *  (((s.z[3]  «  5)  ~  s.z[3])  »  24) 

s.z[4]  =  ((s.z[4]  &  0xffffffffff800000ull)  «  8)  *  (((s.z[4]  «  3)  *  s.z[4])  »  33) 

return  (  s.z[0]  ~  s.z[l]  "  s.z[2]  "  s.z[3]  *  s.z[4]  )  ; 

} 

void  rng : : init(  )  { 

//  the  first  rng  state  s  is  seeded  as  indicated 
//  for  subsequent  rngs  compute  state  s 
//  jumped  ahead  from  the  reference  state 
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//  set  reference  state  for  next  jump  to  zero  state  of  current  rng 


if  (  ref_counter  ==  0  )  { 
seedgenC  ) ; 

jpunpackC  Jump_P[0]  ); 

} 

else 
jumpC  ); 
ref_state  =  s  ; 
id  =  ref_counter++; 

} 

void  rng :: seedgenC  )  { 

//  seed  the  system 

//  set  seed  for  first  generator  (algorithm  specific  to  rng) 

//  the  rng  must  be  in  state  zero,  not  the  seed  state 

s.z[0]  =  seed  *  69069; 

if  (  s.z[0]  <  2  )  s.z[0]  +=  2u; 

s.z[l]  =  s.z[0]  *  69069; 

if  (  s.z[l]  <  8  )  s.z[l]  +=  8u; 

s.z[2]  =  s . z [ 1]  *  69069; 

if  (  s.z[2]  <  16  )  s.z[2]  +=  16u; 

s.z[3]  =  s.z[2]  *  69069; 

if  (  s.z[3]  <  128  )  s.z[3]  +=  128u; 

s.z[4]  =  s.z[3]  *  69069; 

if  (  s.z[4]  <  2048  )  s.z[4]  +=  2048u; 

gen() ; 

} 

void  rng :: jpunpackC  const  unsigned  long  long*  const  Jump_P  )  { 
for  (  int  j  =  0  ;  j  <  5  ;  j++  )  { 

JP[j]  |=  unsignedC  Jump_P[j]  »  32  &  Oxffffffff  ); 

JP[j]  «=  32; 

J P [ j ]  |=  unsignedC  Jump_P[j]  &  Oxffffffff  )  ; 

} 

} 

void  rng:: jumpC  )  { 

//  set  rng  state  s  to  reference  state 

//  compute  temp_state  =  J[0]*s[O]  +  ...  +  J [63] *s [63] 

//  set  s  to  temp_state 

//  now  s  is  jumped  ahead  from  reference  state  by  J 
int  i  ,  j ; 
state  temp_state; 
s  =  ref_state; 

memset(  temp_state.z  ,  0  ,  sizeof (temp_state . z)  ); 
for  (  i  =  0  ;  i  <  64  ;  i++  )  { 
for  (  j  =  0  ;  j  <  5  ;  j++  ) 
if  (  JP[j][i]  ) 

temp_state.z[j]  "=  s . z [ j ] ; 
gen() ; 

} 

s  =  temp_state  ; 

} 

//  number  of  rngs  instantiated 
unsigned  rng: :ref_counter=0; 

//  reference  state  for  jump  computation  is  the 
//  zero  state  of  the  previous  instance 
rng:: state  rng: :ref_state; 
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//  seed  for  first  rng  instance 
unsigned  long  long  rng : : seed=0x!2340fG0 ; 


//  lump  polynomial 

bitset<64>  rng::JP[5]  =  {  0,  0,  0,  0,  0  }; 


//  Tausll3  jump  coefficients 

unsigned  long  long  rng : : lump_P [1] [5]  =  { 

{ 

//  2*2G 

Gx7fd7f8ffc4Gb6103ull  ,  GxGG08a943bd59e7bfull  ,  0xG00bOe2fOe8a5511ull  , 
GxQGGQ50f342c7f5ddull  ,  GxGG0GG0d9bd741 142ull 

//  2*8G 

//  Gx7G2fffbefc7e2115ull  ,  GxGG7fG672769e6G0cull  ,  0xG00aa7dad2018eedull  , 
//  GxGG0G598fll394622ull  ,  GxGG0GGOGG8GGQlGO2ull 

//  2*128 

//  Gx0G0QQ00Q00GO0G10ull  ,  GxGG57eb3124facG97ull  ,  0xG00e88c757dG5b63ull  , 
//  GxOG0O5f36f89389baull  ,  GxGGOGGOG lGGGOGGOGull 

} 

}; 
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C.6  Mersenne  Twister  mt32  RNG  class  rng_jumpjnt32  .h 


//  -*-C++-*- 

//  rng_jump_mt32 .h 

//  Mersenne  Twister  Random  Number  Generator 
//  with  statistically  independent  instances 
#include  <cstring> 

#include  <bitset> 

class  rng  { 

static  const  unsigned  N  =  624; 

struct  state  {  //  state  information 

unsigned  n; 
unsigned  y; 
unsigned  z [N] ; 

}  I 

public: 

//  the  first  instance  is  initialized  from  the  seed 
//  subsequent  instances  are  jumped  ahead  giving  independent  rngs 
rng(  );  //  constructor  (called  with  no  arguments) 

rng(  const  unsigned&  seed  );  //  constructor  (called  with  seed  argument) 

unsigned  gen(  );  //  rng  returns  a  32-bit  unsigned  integer 

static  void  setseed(  const  unsigned&  seed  ) ;  //  call  once  for  all  rngs 
static  const  int  wid  =8;  //  integer  format  width 

private : 

state  s  ;  //  rng  state 

unsigned  id  ;  //  rng  sequential  id  per  instance 

void  init(  );  //  initialize  a  rng 

void  seedgen(  ) ;  //  seed  the  first  rng 

void  jump(  );  //  jump  to  new  state 

void  jpunpack(  const  unsigned*  const  Jump_P  ); 

static  unsigned  seed;  //  seed  for  first  generator 

static  unsigned  ref_counter;  //  number  of  instances 

static  state  ref_state;  //  zero  state  of  newest  instance 

static  unsigned  Jump_P[][N];  //  jump  polynomial  coefficients 

static  const  unsigned  M  =  397; 

static  const  unsigned  MATRIX_A  =  0x9908b0df;  //  constant  vector  a 
static  const  unsigned  UPPER_MASK  =  0x80000000;  //  most  significant  w-r  bits 
static  const  unsigned  LOWER_MASK  =  0x7fffffff;  //  least  significant  r  bits 
static  const  unsigned  CPdeg  =  19937;  //  characteristic  polynomial  degree 
static  const  unsigned  Jdmax  =  CPdeg-1;  //  max  jump  polynomial  degree 
static  bitset<CPdeg>  JP;  //  jump  polynomial 

}; 

rng:: rng (  )  {  //  constructor  (called  with  no  arguments) 

init(  ); 

} 

rng:: rng (  const  unsigned&  seed  )  {  //  constructor  (called  with  seed  argument) 
setseed(  seed  ) ; 
init(  ); 

} 

void  rng : : setseed(  const  unsigned&  seed  )  { 
if  (  ref_counter  ==  0  ) 
rng : : seed  =  seed; 

} 

unsigned  rng::gen(  )  {  //  rng  returns  a  32-bit  unsigned  integer 

unsigned  y; 

static  unsigned  magOl [2] ={0x0 ,  MATRIX_A} ; 

if  (s.n  >=  N)  {  //  generate  N  words  at  one  time 
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int  kk; 

for  (kk=0;kk<N-M;kk++)  { 

s.y  =  (s.z[kk]&UPPER_MASK) | (s.z[kk+l]&LOWER_MASK) ; 
s.z[kk]  =  s.z[kk+M]  “  (s.y  »  1)  '  mag®l[s.y  &  0x1]; 

} 

for  ( ; kk<N-l ;kk++)  { 

s.y  =  (s.z[kk]&UPPER_MASK) | (s.z[kk+l]&LOWER_MASK) ; 
s.z[kk]  =  s.z[kk+(M-N)]  “  (s.y  »  1)  "  mag01[s.y  &  0x1]; 

} 

s.y  =  (s.z[N-l]&UPPER_MASK) | (s.z[0]&LOWER_MASK) ; 
s.z[N-l]  =  s.z[M-l]  '  (s.y  »  1)  "  mag01[s.y  &  0x1]; 
s.n  =  0; 

} 

y  =  s.y  =  s.z[s.n++] ; 

//  Tempering 
y  ~=  (y  »  11); 
y  ~=  (y  «  7)  &  0x9d2c5680; 
y  '=  (y  «  15)  &  0xefc6@0©0; 
y  ~=  (y  »  18) ; 
return  y; 

} 

void  rng : : init(  )  { 

//  the  first  rng  state  s  is  seeded  as  indicated 
//  for  subsequent  rngs  compute  state  s 
//  jumped  ahead  from  the  reference  state 

//  set  reference  state  for  next  jump  to  zero  state  of  current  rng 
if  (  ref_counter  ==  0  )  { 
seedgenC  ) ; 

jpunpackC  Jump_P[0]  ); 

} 

else 

jump(  ); 
ref_state  =  s  ; 
id  =  ref_counter++; 


void  rng :: seedgenC  )  { 

//  seed  the  system 

//  set  seed  for  first  generator  (algorithm  specific  to  rng) 

//  the  rng  must  be  in  state  zero,  not  the  seed  state 
s.z[0]=  seed  &  ©xffffffff; 
for  (s.n=l;  s.n<N;  s.n++)  { 

s.z[s.n]  =  (1812433253  *  (s.z[s.n-l]  "  (s.z[s.n-l]  »  30))  +  s.n); 
s.z[s.n]  &=  Oxffffffff; 

} 

gen() ;  //  for  all  rngs 


void  rng :: jpunpackC  const  unsigned*  const  Jump_P  )  { 
for  (  int  j=N-l;  j;  j —  ) 

(  JP  |=  Jump_P[j]  )  «=  32; 

JP  |=  Jump_P[Q] ; 

} 

void  rng::jump(  )  { 

//  compute  jump  state  t.z[0  ..  (N- 1) ] 

//  from  jump  polynomial  J  and  base  states  sy=s.y 

// 

//  t.z[Q  ]  =  J[0]*sy[©  ]  +  ..  +  J[d]*sy[d  ]  +  ..  +  J  [Jdmax]*sy  [Jdmax  ] 
//  t.z[i  ]  =  J[0]*sy[i  ]  +  ..  +  J[d]*sy[d+i  ]  +  ..  +  J  [Jdmax]  *sy[Jdmax+i  ] 
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//  t.z[N-l]  =  J[0]*sy[N-l]  +  ..  +  J [d] *sy [d+N-1]  +  ..  +  J [Jdmax] *sy[ Jdmax+N-1] 
state  t; 

unsigned  sy [Jdmax+N] ;  //  Jdmax+N  states  0  ..  (Jdmax+N-1)  for  jump 
s  =  ref_state;  //  start  s  at  ref  state 

for  (  int  i=Q  ;  i< Jdmax+N  ;  i++  )  { 

sy[i]  =  s.y  ;  //  save  states  sy[i  in  0  . .  (Jdmax+N-1)] 

gen()  ; 

} 

memset(  t.z  ,  0  ,  sizeof(t.z)  ); 

//  J  term  power  =  d  in  0  . .  Jdmax 
int  i,  k; 

for  (  int  d=®  ;  d<= Jdmax  ;  d++  ) 
if  (  JP[d]  ==  1  ) 

//  i  in  0  . .  (N-l) 

//  k  in  d  . .  (d+N-1)  for  each  d 
//  k  in  0  . .  (Jdmax+N-1)  overall 
for  (  i=®  ,  k=d  ;  i<N  ;  i++  ,  k++  ) 
t.z[i]  ~=  sy[k]  ;  //  k  =  d  +  i 

s  =  t  ; 
s.y  =  s . z [0] ; 
s.n  =  1; 


//  number  of  rngs  instantiated 
unsigned  rng: :ref_counter=® ; 

//  reference  state  for  jump  computation  is  the 
//  zero  state  of  the  previous  instance 
rng:: state  rng: :ref_state; 

//  seed  for  first  rng  instance 
unsigned  rng : : seed=Qxl234Qf®0 ; 

//  Jump  polynomial 

bitsetcrng : : CPdeg>  rng::JP  =  0; 

//  MT32  jump  coefficients 

unsigned  rng : : Jump_P[l] [N]  =  {  //  jump  polynomial 

{ 
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//  2"2® 

8x6df32afe 

8x89f9®a3b 

8x93e769bc 

®xf®6d9493 

®x7®54523c 

®x®f992®b3 

8x4a£627c4 

8xfbl4df97 

®xaeb337®b 

8x7361dc46 

8xcb8el763 

®xfecc9766 

®xeb2ffe48 

8x2bb6d4d8 

®xe9491785 

®xfcd8ee4Q 

®xl94b®bdc 

®x76d8abce 

8x58135637 

®xaa®2cbac 

8xld558c63 

®x6Qaceflf 

®xcl43f425 

®x4dl7cc6® 

8x7®8e558d 

8xl987el2c 

8xe6e6c8fd 

®xc7857b86 

8x4866c355 

8x45662187 

8x4b828®6a 

8x47f25dl3 

®xl41a6e5c 

8xb2aae41a 

8xe61c276d 

®x3897®89d 

®xd5cc3Q4c 

8x6bed6442 

8x9cc64el® 

®xc2568953 

®x6eac7®44 

®xl3®acc82 

8x993ae794 

®x94a®e6a2 

®x®38ca3®l 

8xc7£af2d9 

®x9a6c4eb6 

8x91ed8767 

®xa6d6b®a9 

8x38£leb3c 

®xbb56ac97 

8x8®f6459® 

8x85flf®b2 

8x4®f7e374 

8xc854c216 

8x8fell4£5 

®xfab56cel 

®x63f8®9ed 

®xd3efa®8c 

®x6ff®8875 

8x8e377962 

®xed®e7c75 

®x4689d®46 

8xf764e416 

®x2256®9®e 

8x82®5a3®b 

8xacf26d97 

®xd97ca®4e 

8x333b4e9f 

8x476b987b 

8xc568e2bf 

®xa2®ee8cd 

8x54435fc7 

®x4d8f8d6e 

®x4eec951® 

®x74118fe6 

®xce5894bb 

®x963913®f 


Qx7eab®cac 
8x6769 l£cd 
8x6®683£55 
SxeccdefSl 
8x42997177 
8x7b9985ca 
8x2698772c 
8x®c5144c2 
8x4fe7a88b 
8x6a357®4® 
8x®2d67f85 
8x5963723d 
8x21f46276 
8x83a48bae 
8xb®9b966£ 
8xf33d987c 
8x9aalf5d® 
8x7518f2f4 
8x27f8726c 
8xc®2ell68 
8x29f®cab4 
8xcc382c9f 
8x®23c3573 
®xeabll®7c 
8x35c52964 
8xe4d54472 
8x89281233 
8xed67e22e 
8x9ec898fa 
8x7b2df£59 
8x2bac®545 
8x42f71ffb 
8x55bfcdd6 
8xbad5bldb 
8xd92d83cd 
Qxbal48634 
8xb818d6ea 
8xc®371b44 
8xd874afal 
8xdeea32ed 
8x73211586 
8xel6adee2 
8xdlda7674 
8x5b877e7e 
8x825bl®la 
8xf9eaa8d4 
8xl£1199e6 
8x2a85b8d7 
8xf®a8e819 
8xa2aec7b£ 
8xd98189f2 
8xbd7c8a75 
8x7c99b97e 
8xc9b6bc5e 
8xdba749dd 
8x8596b44® 
8xdaa587b5 
8x®61bd7c3 
8xc22ce857 
8xc4bl9427 
8x7eefc4a3 
8xc2ded347 
Qx2a7®7427 
8xblfd7da7 
Qxf 1295cQ4 
8xl4b2afa5 
8x783delf9 
8x48e35®39 
8xd9ea551e 
8x62dee751 
Qx636ade37 
8xl3735ce5 
8xdf32c£3a 
8x8e42994d 
®xb®571327 
Sxcddbcada 
8x®151b56c 
8xaabb9e4d 


8xbefa4d71 
8xc612f2 lb 
8xf86e41bc 
8x7e2be6a2 
8xl787d5®9 
8xc6b®c6cl 
8x735adb®8 
8x3fl6b27e 
8x®3d8bfd® 
8xad3442f4 
8xbl9693f8 
8x3cb73c9b 
8xd®b938df 
8xle48al37 
8xf6c96f54 
8xc893f8ae 
8x59b7fbc5 
8x54d41dl6 
8xa8d64246 
8x52722adc 
8x2378648b 
8xd8188b7a 
8x9ec7da44 
8x6483cc81 
8xb®fb49eb 
8x453237dc 
8x62571b®2 
8xb8b7dac8 
8xeca55312 
8x7d3c8®71 
8x8ab681f3 
8x®9dfafd6 
8xd2fcf28a 
8xe5c7e®d9 
8xcfcd327£ 
8xb3d53f2a 
8x835cb®d4 
8x35f273c4 
8xl98d2a62 
8xb8346176 
8xl644b661 
8x®83824e4 
8x85b818fa 
8x333473ad 
8x691d2d45 
8x5eb7e2fd 
8x8278379£ 
8x53cfa372 
8xle®aba84 
8xd2e73d92 
8x5aa6591c 
8x66562d6d 
8x4b®113e3 
8x8fl4a77e 
8x997737a £ 
8x678e7f4f 
8xd8d34788 
8x2d®679e2 
8xflacl25c 
8xe2c®8c3c 
8x7275f6e7 
8x9e3d®dc8 
8xbe623eed 
8xd®24e345 
8xel4f4cl5 
8xdac6984d 
8x®2cc9618 
8x®cda8881 
®x®68clcfe 
8xb8c2e9af 
8xb88eb8a® 
8xee967dd4 
8x4a3d®22c 
8x73bb383c 
8xacealb64 
8x7clla9cd 
8xde8117c£ 
8x8c5658d9 


8xadc4c8b2 

8x58846827 

8x5d567f3c 

8x3c£24cl3 

®xdebc®78a 

8x4dd5282c 

8x2ee6afdc 

8xcl8df4b3 

8xbb7da76e 

8xcf778ab® 

8xc4c7246b 

8xae9efa7b 

8x45657132 

8x64cd3bc7 

8x923a8c®d 

8xl2486eel 

8xl563fdbc 

8x9c6e661d 

8x8d79c®54 

8x4f5a25®5 

8xa6e®63£7 

8xa8b82496 

8xbc47fffe 

8x3d868fl9 

®xfb®59a®2 

8xb68937df 

8xf®8df485 

8x®17fb4e6 

8xd8156db6 

8xc79def33 

®xc®8228bd 

8x6b93afc7 

8xl4ee8ed7 

8xffe7b659 

8xb®83cd35 

8x488bl74c 

8x83226722 

8xdc85559® 

8xl3d92dcf 

8xbdca4be9 

8x64914aa8 

8x2adef6ba 

8x3fbl84c4 

8x4899a4fe 

8x2c79£e2e 

8x8598ba54 

8xc2984f88 

Qx56d®7626 

8xadbb42d8 

8x79a525bd 

8x92485e89 

8xacb28d28 

8x5628bafd 

8x6e8f22ea 

8xl2651c8a 

8xl5866f46 

8xacl66b9b 

Qx58d2476e 

8x97b49fe6 

8xe9bfl881 

8xe9c9f32b 

8x498271fb 

8xfed455£f 

8xedbb9dd2 

8xd®a26a57 

8xe5cel8d8 

8x45984d87 

8xe68481ef 

8x9181el6c 

8x3289e7e2 

8x82ba431a 

8x21d58e76 

8x4c48e3e4 

8x3d7f9118 

8xdb9baa85 

®xecdfl®8a 

8xcl69®4®9 

8x9cla586b 


8x23f269fl 

8xea84d6f9 

8x2elf®91d 

8xlc9f44db 

8x5aaeb386 

8x48359cfl 

8x7a5a6c81 

8x6eddc®cf 

8xd37ba3d4 

8x®82d3bfe 

8xc6ccda®l 

8xf4e84clf 

8xbllecaa4 

8xld79c98d 

8x47elc778 

8x61de23f5 

8xa®b8c39c 

8x®9cf99fc 

8xel35ce55 

8x522b22f9 

8x34e436e7 

8xc81e5237 

8xd37b9edl 

8x5be5a795 

8xc74e®8d9 

8xd7878393 

8x®956b5d9 

8x4b9883e5 

8xa5ba88®7 

8xl61a®9fa 

8xc482d®a5 

8xbe8aeb®l 

8xlc93d772 

8xbae4c34e 

8x3821f57e 

8xe®91d5f3 

8x6796e622 

8xff734de5 

8x938539c® 

8xfe®ba665 

8x38fc5881 

8x66864277 

8x9d447al7 

8xffe®5£c8 

8xf8678882 

8x45856d6d 

8xde7a97ea 

8x662ec41b 

8x7e64ae3a 

8x9f8e4e®f 

8x645321ac 

8xcelfele7 

8x8ba7e®42 

8x4bf2d®14 

8x82245826 

8x3bb4e9f9 

8x24f64548 

8x34826f®4 

8xl2a8187b 

8x9858®b2f 

8xfd2c46db 

8x78ff97ca 

8x53698eac 

8x79c99b3c 

8x®ffcaal2 

8xcbbf27b8 

8xf558eb3f 

8x25921fd® 

8xl®81af6c 

8x719e9b4a 

8x272del81 

8x65a7d®16 

8xl24b242e 

8xc399el68 

8xa327e725 

8x64c83e51 

8x7d282b3c 

8xde395cb® 


®xf®955951 

®xa®24724a 

8x692a7e93 

8x3fb5d8a4 

8x931fd6a6 

Qx32e3cb72 

8x®898faf7 

8xad5a221c 

8x4f®b4778 

8xl2915fe6 

8xle7a413b 

8xf2bb®2aa 

8xcaa88258 

8x28fdalll 

8x53348854 

8x2c58648b 

8x84ceeaba 

8x62967b52 

8x8cal6ff9 

8xac4f438® 

8xcefc457a 

8x83671917 

8x2a9e7eb® 

8x7a3aa33e 

8xa25elc8d 

8x7c88fe2e 

8xa4f3f91c 

8xb57b6432 

8x4dcc4554 

8xdfbflfb2 

8xa76f5824 

8x79bc®dlf 

8xcl29888e 

8xe3a22®59 

8x372ca®22 

8x7589c8bb 

8x2882eb5® 

8x28848dea 

8xce691®92 

8x3392ela4 

8xec3f5ale 

8x®8e7ba58 

8x7f4f38ac 

8x3148b®f2 

Qxdd7fl37e 

8x2849d7a3 

8x6ae6b57a 

Qx3b8b74b7 

8x®49bc863 

8x266d®ca2 

8xa92743e7 

Qx7al38a®l 

8x9ab5d848 

8xab61d46d 

8x367edl9® 

8x558272a2 

8x24bled72 

8x686c948f 

8x387f36®2 

8xc58d42b8 

8x63eb94fb 

8x4aabl51b 

8xb848318d 

8xl3b®e5e5 

8x83b8c4f2 

8x4fb7d®lc 

8x®5ef6841 

8x9c6be6d5 

8xd2d2dd81 

8x53615727 

8xafd623c8 

8x®8b6e7b® 

8xl4dlaade 

8xc6alb279 

8xaa682a2b 

8xe66ccf5f 

8xcd48ee59 

8x38841ec8 


8xal9bc77c 
8xlle®ae46 
8x2359e£12 
8xl44d3faf 
8xa7b8ddfe 
8xa883®761 
8xbaflc4a2 
8x5145bed4 
8x8931dl9® 
8xfla38cbl 
8xbd6®666® 
Qx3aQal227 
8x82dae28£ 
8xe8f98e31 
8xe3f7fb74 
8x7a58d21£ 
8x2484537® 
8xfee588bl 
8x661a88e£ 
8xb974fle9 
8xa7187b46 
8xbl7db7bd 
8xaddceca4 
8xe239e614 
8x3e33clle 
8xd23d2a®b 
8x46a631e£ 
8x38f37aa6 
Qx799e57d3 
8x914d5flb 
8xf52bb5f3 
8x4e843718 
8x61885675 
8xa51f6£a7 
8xbd®83b47 
8xl5cedae® 
8xc8b4968c 
8x5f9df6c7 
8xbf2dced7 
8xb6d7ee®e 
8x33782633 
8xf4bfllb9 
Qx33aad744 
8x®3f4d5bc 
8x648ec®2b 
®x®458e®3b 
8xe688ce9d 
8x®dl8bea3 
8xbebfale4 
8xb3f8b9fb 
8x811d2e95 
8xd63ea468 
8xe3e6b974 
8xd®d4394® 
8x698468d4 
8xlcb36315 
8xd®daff29 
8x®99e8172 
8x554673ba 
8x21c93ca5 
8x5873a894 
8x7218d914 
8xf3fcb7b® 
Qx51ca5dd3 
8xe4887c96 
8xd53bc9b8 
8x88belee7 
8x38b79457 
8x3bc888bb 
8xa85aa982 
8xlac6475a 
8x631e38d® 
8xffa63459 
8x377 lba®4 
8x8efc4ad8 
8x6353efd9 
8xd68®2857 
8x43368adc 


8xccf®ld2f 
8x7flc4ca3 
8x8b58b2c® 
8x328892e9 
8xe®3e657b 
8x5c433f85 
8xd62a8fea 
8xad742fe2 
8xc4c93cbd 
8xf®69cd8® 
8xal8a553a 
8xeefe6b5a 
8xf7de898® 
8xc389e322 
8x518245ca 
8x913b6d92 
8x4c78236c 
8x2a734b9f 
8xlaedl78c 
8xc8e41ae8 
8x2856d593 
8xa82f8b®2 
8x91cad®71 
8x8c84cff5 
8x4baff421 
8x98fb33ec 
8x7cbe51el 
8xf4£c73d3 
8xd3413725 
8x8682 lf42 
8xe®£55c25 
8x5c8cl3£c 
8x27e22987 
8x®5816c33 
8x®dclaa8d 
8xc42d88e6 
8x6dc5aab8 
8xfa4elcd5 
®xf®43abdd 
8xe7287f2e 
8x8f3ccff3 
8x82a9435b 
8xb744aa59 
8xb7166e53 
8x2458138c 
8xcc72979® 
8xa869aa4c 
8x2c2fl9f® 
8x8cll®c39 
8xl2cbl®8® 
8xl9dfld2b 
8xlb4497a5 
8x77932da8 
8xe28d®74b 
8xf4bb2137 
8xc9baa847 
8x4a52eall 
8x41982863 
®x®2442fae 
8x22ac3484 
8x8758e968 
8x4d818e68 
8x48dc3283 
8xdae497e3 
8xa7a6a41b 
8x®2d9ce36 
8x39a2783f 
8x54994e56 
8xb3fef358 
8xe79f®c32 
8x2c6beaee 
8x8dl639a9 
8x4f3c®862 
8xd533£333 
8xcae6bdc3 
8xf98d6ccd 
8xfl2318f3 
8x88888881 
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//  2~imm 

//  ®xc®4347b3  ,  8x568188e5  ,  8xf6cld4d®  ,  8xdd9ff6f9  ,  ®xle®7ec3a  ,  8x27c9cd71  ,  8xdb5f3e52  ,  ®x8fd2c779 

//  ®xf33dlf®f  ,  8x35®a5ea5  ,  8x4b®98eb2  ,  Sx9c2e5c2f  ,  Qx2cfl9af8  ,  ®xf2515f4®  ,  8x8fb9cfa®  ,  ®xfaa5a36c 

//  ®x7a47cd4a  ,  8xb22e9739  ,  ®x8e34Q824  ,  8xb515d2b£  ,  8x716e31®3  ,  8xda®a315f  ,  ®xfcclcce4  ,  8x4d2®c264 

//  8xf9d6bc5b  ,  ®xea®®f6®c  ,  ®xdfdl3eee  ,  8x94a9d212  ,  Qx85fb2682  ,  8xe8al4ffd  ,  ®x816®533c  ,  8x512c3ec3 

//  8x48df671c  ,  8xcd5b5a7e  ,  8xd76a8f57  ,  8xeb84a72b  ,  8x3e71£91®  ,  8xe44f819®  ,  ®x5521bb®c  ,  8xb7bb65bl 

//  8x7c8a942e  ,  ®xbd7®4511  ,  8x5el5e781  ,  ®xl2e3ccle  ,  Qxa736a®e3  ,  8x437998a9  ,  ®xaf362®81  ,  8x42a8e62b 

//  Qx®833ad3e  ,  ®x54dedalc  ,  8x4f®5edl3  ,  8x2a®a7124  ,  Qxal9449ba  ,  8x45®3®aad  ,  8x5c84d8d7  ,  ®xcf6®fdla 

//  ®x891db7cf  ,  8x96c26ee3  ,  8xaca98ac8  ,  8xa71£a59e  ,  Qxe631c7b3  ,  ®x387elb29  ,  8x2e®2a376  ,  8xb98cc62c 

//  8x79®624f4  ,  8xcle64a2e  ,  ®xfl66®529  ,  8x595e5fbe  ,  8xc334el47  ,  Qx4®2eb5bl  ,  8x7ceacc2f  ,  8x3666b776 

//  8xa3dlf2c5  ,  8x22683362  ,  8x71361ed6  ,  8xaa2df833  ,  8x8636a3b6  ,  8x7c711139  ,  8xee71425b  ,  8x631ed752 

//  8x7c2131®9  ,  8x937e®18c  ,  8xll2f65ec  ,  8x2b94197a  ,  8xcb6944b6  ,  8xb2ea7®12  ,  8xcaa®edd4  ,  8xlaa56252 

//  8x2abal869  ,  8x852a69c5  ,  8x278f8cda  ,  8x4519f6ee  ,  ®x®de®3259  ,  8x66687b7e  ,  8xa2b6fa7b  ,  8xal74888e 

//  8xe4e4daa4  ,  8x71cb36®4  ,  8x87987d88  ,  8xldld8543  ,  8x61fb£37c  ,  8x72a7bbda  ,  8xdf38fa82  ,  8x48e9ac4® 

//  8x67e35fed  ,  8x823676c®  ,  8x3c3af93f  ,  8xd793b469  ,  8xf89dl2ad  ,  8x48c6a78£  ,  8x28783194  ,  8xaadc641c 

//  8xc8bl63f5  ,  8x7a£249®9  ,  8xfbba7a®8  ,  8x91a76daa  ,  8x8118928b  ,  8x2a72596a  ,  8xe®c7678b  ,  8xe25bbb®3 

//  8x9el38f3£  ,  8x286de6e2  ,  8xa37cc871  ,  8xde3f8c®8  ,  8xcdd®8cc2  ,  8x965c61aa  ,  8xfd2dl6a3  ,  8x52aebd6b 

//  8x59865fc®  ,  8xfb7d5c5c  ,  8xd27b57£e  ,  8xl9d345d6  ,  8x8c255798  ,  8xa62a62ba  ,  8xec3e7b3b  ,  8x85b5637d 

//  8x214ae6eb  ,  8x®9186b86  ,  8xl5c513£5  ,  8xbc526b75  ,  8xl37dllb2  ,  8x22e2c7f4  ,  8xf8c3819f  ,  8xc34a821e 

//  8xcf24875a  ,  8x242efcf®  ,  8x5892b749  ,  ®xccc®8afc  ,  8x5852bf3d  ,  8x®bl4836b  ,  8x88ad786d  ,  8x45a76d4e 

//  8x72fdbcl4  ,  8xe549d211  ,  8xe56c3d7e  ,  8xf5f7d5®7  ,  8xdlface87  ,  8x6668d5cl  ,  8xbe44a38e  ,  8x52847163 

//  8x®d6ebcb®  ,  8xa2395417  ,  8x®fb22d85  ,  8xdll94dl4  ,  8x6dcaf675  ,  8x6953eec4  ,  8x86313d46  ,  8xbb291c7a 

//  8x56b57c64  ,  8x8d243f64  ,  8x2a8f985e  ,  8x6522af6b  ,  8xbf6ce5ee  ,  8xl45387f3  ,  8xfa7eecf7  ,  8xee544df4 

//  8xble®37a2  ,  8xaflla®d5  ,  8xd69d45b5  ,  8xed78f57c  ,  8x854e244a  ,  8xb4e227®8  ,  ®xc®49eae®  ,  8xlaa518d5 

//  8xce98c52£  ,  8x512d47df  ,  8xeb2e2776  ,  8xbd273cba  ,  8xa6a®7e7f  ,  8x7bf®5af2  ,  8x81625574  ,  8x78c859c£ 

//  8x6d5a36f8  ,  8xala2bb8a  ,  8xf83a483a  ,  8x2886f2bl  ,  8xbf3a56e®  ,  8x55df7b2®  ,  8xe6872®31  ,  8xace8583e 

//  8xef93e5®5  ,  8xc65a7bfe  ,  8x®78el2f2  ,  8xd6ed5d2b  ,  8xcfclaa®f  ,  8x5341a38d  ,  8x37be5beb  ,  8x7d3d8871 

//  8x4c5f4ffa  ,  8xa5962cl9  ,  8xecfda25d  ,  8xb4c9d71e  ,  8x®cbcf349  ,  8x56f8812b  ,  8x4f46352c  ,  8xc673a535 

//  8xa6efa2a2  ,  8xe8c5e384  ,  8x48flee45  ,  8x7£d39827  ,  8xlc5cl®ld  ,  8xd®6818dl  ,  8xffa5f89a  ,  8x9353e4ce 

//  8x39b5c7f2  ,  8x93e6b28c  ,  8xf32818c2  ,  8x418d424a  ,  8xcf33e3db  ,  ®xfld93ffc  ,  8x2371abd2  ,  8xe6794e9d 

//  8x388db2fa  ,  8x594a8886  ,  8x82367b5d  ,  8xe6552e97  ,  8x3fdclea®  ,  8xf385357d  ,  8x7e31868a  ,  8x5d27cc®5 

//  8x32c69f31  ,  8x8167f497  ,  8x7e6al761  ,  8x88al3229  ,  8xb776533a  ,  8x44949f6e  ,  8x2a5bb57e  ,  8xel32d87b 

//  8x8ed82aff  ,  8xac3366f6  ,  8x2fc®d961  ,  8xdc8ca748  ,  8xc963d571  ,  8xcdl®d6c7  ,  8x4af®d83a  ,  8x22ad425e 

//  8x9438338a  ,  8x572 54ab7  ,  8xed686c38  ,  8xb262df3d  ,  8xe9b2adcb  ,  8xb94d5889  ,  ®xe®6d®bdl  ,  8x2e488dab 

//  8xb®al52cd  ,  8x367a3f21  ,  8xf5cafd56  ,  8x9b7c8c®a  ,  8xba4c3616  ,  8xf8e64f®8  ,  8x68c84233  ,  8x9d677ac6 

//  8x96692 18f  ,  8xdbf®68e4  ,  8xb585eb56  ,  8xl£cll3a7  ,  8x437cd7d8  ,  8xlc4df6b®  ,  8x5ac7d45d  ,  8xae884cab 

//  8x53874ef5  ,  8x3e673518  ,  8x9919el41  ,  8xf56de®48  ,  8x®6556d48  ,  8x38aedadd  ,  8xaf72a9d4  ,  8xb7d8fbl5 

//  8xflcl5b6d  ,  8x33blc23b  ,  8x4fc4369®  ,  8xf5ecl2ef  ,  8xl7ef5®42  ,  8x58cf794£  ,  8x9f354ee9  ,  8xb334e323 

//  8xcc8695f8  ,  8x3dd367ba  ,  8x6826el89  ,  8x8bb7cf81  ,  8x742db792  ,  8xa7a428ba  ,  8xb281c923  ,  8x®c2fe£f£ 

//  8x989e7dc2  ,  8xef68de3d  ,  8xc®a63dd8  ,  8xfflb2a44  ,  8xf235c381  ,  8xel66e49e  ,  8xb38176ba  ,  8x2872d89e 

//  8x418c2®9®  ,  8x5e89e®d9  ,  8xaff5615®  ,  8xe9c32afa  ,  ®xl®513bae  ,  8x2dl3faa®  ,  8xl25511f8  ,  8x9b7d2fl6 

//  8xf6151b®3  ,  8x45d4e4ab  ,  8x75926964  ,  8x24a5e®19  ,  8x9eb8ef£2  ,  8x51cfb®5d  ,  8xa682eb4b  ,  8xf7ad94a6 

//  8x7732df95  ,  8xbd3dbb29  ,  8xd378b2e2  ,  8xc85828d2  ,  8x73c®d682  ,  8x25d61874  ,  8x78bel45b  ,  8x4ed37f82 

//  8x4dbfecb5  ,  8xl4c8148f  ,  8xd®3fd296  ,  8xef294e91  ,  8x78c44a®f  ,  8xa438e693  ,  8x®8ala573  ,  8x414bcbec 

//  8x8121e35b  ,  8x582585c3  ,  8x3895abb4  ,  8xb71d788e  ,  8xe®271bf2  ,  8xa95ae31b  ,  8xcdafd26f  ,  ®xf582affl 

//  8x18896433  ,  8x®la23922  ,  8xddd4233®  ,  8xaab93eda  ,  8xe6f98289  ,  8xb642ed7b  ,  8xlb394182  ,  8xdd91cad® 

//  8x3fb9b8e8  ,  8x881af9e2  ,  8x8d9a92bb  ,  8x716acldf  ,  8xf99489f®  ,  8x32c58®75  ,  8xedl4da69  ,  8x9f2dlb8® 

//  8x45de2ebc  ,  8x96ab3b86  ,  8xddfdf975  ,  8xc877a977  ,  8x®8b526fa  ,  8x255baafe  ,  8x6bffc®72  ,  8xbb7d639® 

//  8x5f25a865  ,  8xb82d®7b5  ,  8xe73dbe9e  ,  8x8cl785a8  ,  8x7ell957e  ,  8x9b538b®4  ,  8xf8f46721  ,  8xclbd6be3 

//  8x98f6el38  ,  8x365eclcd  ,  8x4c®c8ad8  ,  8xa9elb736  ,  8xcc77abbb  ,  8xa43965b4  ,  8x5ale88c7  ,  8x3798ebb3 

//  8x88aalfe3  ,  8xdf33b33a  ,  8x859c42bb  ,  8x71b6da3c  ,  8xec36ce2c  ,  8x27a428®9  ,  8x49bc6a5b  ,  8x9f®8adbd 

//  8xl833ba9a  ,  8x5c®3alc7  ,  8x88f86aa9  ,  8xd5413e5®  ,  8x3439cde®  ,  8x2ec57fb®  ,  8x7913997d  ,  8xc2alfll7 

//  8xl813c975  ,  8x53694a®5  ,  ®xcfd®8eld  ,  8xe66fl93b  ,  8x513254£f  ,  8xf5f6f®b8  ,  8xb23c7f98  ,  8xbdc389fl 

//  8xaf3cb7c2  ,  8x®3e469a2  ,  8x®961a615  ,  8xla862e43  ,  8xb934bllc  ,  8x5b88247b  ,  8x286b43f8  ,  8x4e®bl3de 

//  8x272b8f97  ,  8xa424fd£8  ,  8x4367cc98  ,  8xc2f6382a  ,  8x3c®a52ef  ,  8xc28827f8  ,  8x985c6acb  ,  8x443d4ead 

//  8xf®dd4d6d  ,  8xl6affa44  ,  8x63bb3b41  ,  8x71311b42  ,  8x4b®f34d3  ,  8xa389831d  ,  8x9d®24dcf  ,  8x9362a49d 

//  8x7bd2513e  ,  8x218a8f74  ,  8xe9daac®8  ,  8xae9a9a4®  ,  8x85e838eb  ,  8xll5fb72£  ,  8xed69c4da  ,  8x®193c285 

//  ®xdab®6bed  ,  8x4a898ba®  ,  8xede875c®  ,  8x4686488a  ,  8x58157a£l  ,  8x4fd5a859  ,  8xca7f36ce  ,  8x85a371c7 

//  8x91798f 14  ,  8x58c®b61b  ,  8xb5c2d824  ,  8x94883b9a  ,  8xb91bb396  ,  8xc55fabf2  ,  8x6fl8c384  ,  8x®lbccb78 

//  8xla533bb3  ,  8xl38633ef  ,  8x91dl424e  ,  8x84f8788b  ,  8xb9cdd®5b  ,  8xlb7e4755  ,  8xb687c92e  ,  8xlfbf7fd® 

//  8xel79®294  ,  8xdbdc9182  ,  8x4e5bcl39  ,  8x8a7237a3  ,  8xfcdbd22a  ,  8x818add97  ,  8xf7654b2b  ,  8xe537c88c 

//  8x9dce8147  ,  8x8fflfaf®  ,  8x7dbc4e95  ,  8x23e®6c41  ,  8xcb9855e6  ,  8xea861e33  ,  8x68c6987a  ,  8xb9d395c6 

//  8xfel®9ca3  ,  8x365be®e6  ,  8x35cba7e3  ,  8xf6f88bf9  ,  8xc5484814  ,  8x937bfe73  ,  8x®9423d61  ,  8x9fc35ae8 

//  8x2bfdc9d4  ,  8xll5485b4  ,  8xfd33d922  ,  8xca28a68e  ,  8x76556bc®  ,  8xlebla222  ,  8x3c6692ef  ,  8xe2e93681 

//  8xbc®88f9b  ,  8xa5843c62  ,  8xf53831ce  ,  8x2d7el2b£  ,  8xla®e3dd7  ,  8x3e®d758c  ,  8xe8b9a962  ,  8x8a6clfcd 

//  8x7a685989  ,  8x29289c9®  ,  8x9dbc3f79  ,  8x8e3c7558  ,  8x38695118  ,  8xlc97111e  ,  8xad2e2bfd  ,  8xl94ac45c 

//  8xce869454  ,  8xlel54de8  ,  8xl788db54  ,  8x6fa31781  ,  8xb55a5fl3  ,  8xa2adecbe  ,  8xe4789c72  ,  8x748d®8b8 

//  8xdabe3974  ,  8xb3a81a52  ,  8x2f7683e4  ,  8xb6e583c5  ,  8x5674d357  ,  8x43786b96  ,  8xe2ceaf61  ,  8x48fa5e8a 

//  8xaa733f96  ,  8x7266bdaa  ,  8x526c6988  ,  8xe3174696  ,  8x3fl®f22b  ,  8x6f5f33bl  ,  8x®d8ebcf9  ,  ®x®9554c®£ 

//  ®xe2499364  ,  8x2cd®e56d  ,  8xe4748882  ,  8xe38badc2  ,  8x9f4541d2  ,  8xfee81859  ,  8xf7dabd8®  ,  8x64cc38b7 

//  8xflb97277  ,  8x98911a43  ,  8x822f651c  ,  8xll6f61d6  ,  8x2fba7e44  ,  8x®ca56638  ,  8x®al398f6  ,  8xdf49alll 

//  8xbble3d72  ,  8xfe269c2e  ,  8xbf48b795  ,  8x7692c6c9  ,  8x7fee9bc®  ,  8xa449869a  ,  8xfef9e8a4  ,  8xb52abc3b 

//  8xd6bldf23  ,  8x87bc8b43  ,  8x7cae9ca7  ,  8xalfa251®  ,  8xed25edc4  ,  8x719db881  ,  8x3f92f8dl  ,  8x84f7ff7a 

//  8xd34619b2  ,  8x87889643  ,  8x91238e94  ,  8x3bbl3b8®  ,  8x®2a78c9b  ,  8xc672546a  ,  8x5e2d74®5  ,  8x85fdc®fl 

//  ®x®585438e  ,  8x7925ea85  ,  8x®4f79a64  ,  8xe6d7f86b  ,  8xb3cf44a7  ,  8xc451a®5a  ,  8x9e4cl®fd  ,  8xl828836e 

//  ®x®2325e®2  ,  8x2e3981d9  ,  8xe2857alc  ,  8x73f48aal  ,  8xd8b49a2e  ,  8xld4b68cd  ,  8x846ee8b3  ,  8x58e®efea 

//  8x98dfd®93  ,  8xb9457b66  ,  8x4b3489bc  ,  8xl48fc279  ,  8xd®ebf34d  ,  8xc®63a952  ,  8xecd57eae  ,  8x78e5bda4 

//  8xcb2d88ba  ,  8x2bc5383f  ,  8xf88debe3  ,  8x88e614df  ,  8xf729eb48  ,  8xf792b83c  ,  8x8e2796b4  ,  8x33e4b6eb 

//  8xbfcfl784  ,  8x2439dc92  ,  8x3ebabbbd  ,  8x79143883  ,  8x452321e®  ,  8x495ale4b  ,  8x5dd7876a  ,  8x88888881 

} 

}; 
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C.7  Mersenne  Twister  mt64  RNG  class  rng_jumpjnt64  .h 


//  -*-C++-*- 

//  rng_jump_mt64 .h 

//  Mersenne  Twister  Random  Number  Generator 
//  with  statistically  independent  instances 
#include  <cstring> 

#include  <bitset> 

class  rng  { 

static  const  unsigned  N  =  312; 

struct  state  {  //  state  information 

unsigned  n; 
unsigned  long  long  y; 
unsigned  long  long  z [N] ; 

}  I 

public: 

//  the  first  instance  is  initialized  from  the  seed 
//  subsequent  instances  are  jumped  ahead  giving  independent  rngs 
rng(  );  //  constructor  (called  with  no  arguments) 

rng(  const  unsigned  long  long&  seed  );  //  constructor  (called  with  seed  argument) 
unsigned  long  long  gen(  );  //  rng  returns  a  32-bit  unsigned  integer 

static  void  setseed(  const  unsigned  long  long&  seed  ) ;  //  call  once  for  all  rngs 
static  const  int  wid  =16;  //  integer  format  width 

private : 

state  s  ;  //  rng  state 

unsigned  id  ;  //  rng  sequential  id  per  instance 

void  init(  );  //  initialize  a  rng 

void  seedgen(  ) ;  //  seed  the  first  rng 

void  jump(  );  //  jump  to  new  state 

void  jpunpack(  const  unsigned  long  long*  const  Jump_P  ); 

static  unsigned  long  long  seed;  //  seed  for  first  generator 

static  unsigned  ref_counter;  //  number  of  instances 

static  state  ref_state;  //  zero  state  of  newest  instance 

static  unsigned  long  long  Jump_P[][N];  //  jump  polynomial  coefficients 

static  const  unsigned  int  NN  =  312;  //  rng  parameters 

static  const  unsigned  int  MM  =  156; 

static  const  unsigned  long  long  int  MATRIX_A  =  0xB5026F5AA96619E9ull ; 

static  const  unsigned  long  long  int  UM  =  0XFFFFFFFF80000000U11 ; 

static  const  unsigned  long  long  int  LM  =  0x0QQ00Q0Q7FFFFFFFull ; 

static  const  unsigned  CPdeg  =  19937;  //  characteristic  polynomial  degree 
static  const  unsigned  Jdmax  =  CPdeg-1;  //  max  jump  polynomial  degree 
static  bitset<CPdeg>  JP;  //  jump  polynomial 

}; 

rng:: rng (  )  {  //  constructor  (called  with  no  arguments) 

init(  ); 

} 

rng:: rng (  const  unsigned  long  long&  seed  )  {  //  constructor  (called  with  seed  argument) 
setseed(  seed  ) ; 
init(  ); 

} 

void  rng : : setseed(  const  unsigned  long  long&  seed  )  { 
if  (  ref_counter  ==  0  ) 
rng : : seed  =  seed; 

} 

unsigned  long  long  rng::gen(  )  {  //  rng  returns  a  32-bit  unsigned  integer 

unsigned  long  long  int  y; 

static  unsigned  long  long  int  magQl [2]={0xQull ,  MATRIX_A}; 
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//  generate  N  words  at  one  time 


if  (s.n  >=  NN) 

{ 

int  kk; 

for  (  kk  =  0  ;  kk  <  NN  -  MM  ;  kk++  )  { 

s.y  =  (  s.z[kk]  &  UM  )  |  (  s.z[kk+l]  &  LM  ) ; 

s.z[kk]  =  s.z[kk+MM]  "  (s.y  »  1)  ~  mag01[(int)  (s.y  &  Oxlull)]; 

} 

for  (  ;  kk  <  NN  -  1  ;  kk++  )  { 

s.y  =  (  s.z[kk]  &  UM  )  |  (  s.z[kk+l]  &  LM  ) ; 

s.z[kk]  =  s.z[kk+(MM-NN)]  "  (s.y  »  1)  "  mag01[(int)  (s.y  &  Oxlull)]; 

} 

s.y  =  (  s.z[NN-l]  &  UM  )  |  (  s.z[0]  &  LM  ) ; 

s.z[NN-l]  =  s.z[MM-l]  "  (s.y  »  1)  ~  mag01[  (int)  (s.y  &  Oxlull)]; 
s.n  =  0; 

} 

y  =  s.y  =  s.z[s.n++] ; 

//  Tempering 

y  (y  »  29)  &  Ox5555555555555555ull ; 
y  ~=  (y  «  17)  &  0x71d67fffeda60000ull ; 
y  ~=  (y  «  37)  &  0xfff7eee000000000ull ; 
y  *=  (y  »  43)  ; 
return  y; 


void  rng : : init(  )  { 

//  the  first  rng  state  s  is  seeded  as  indicated 
//  for  subsequent  rngs  compute  state  s 
//  jumped  ahead  from  the  reference  state 

//  set  reference  state  for  next  jump  to  zero  state  of  current  rng 
if  (  ref_counter  ==  0  )  { 
seedgen(  ) ; 

jpunpack(  Jump_P[0]  ); 

} 

else 
jump(  ); 
ref_state  =  s  ; 
id  =  ref_counter++; 

} 

void  rng : : seedgen(  )  { 

//  seed  the  system 

//  set  seed  for  first  generator  (algorithm  specific  to  rng) 

//  the  rng  must  be  in  state  zero,  not  the  seed  state 
s.z[Q]=  seed  ; 

for  (  s.n  =  1  ;  s.n<NN  ;  s.n++) 

s.z[s.n]  =  (6364136223846793005ull  *  (s.z[s.n-l]  ~  (s.z[s.n-l]  »  62))  +  s.n); 
gen() ; 

} 

void  rng : : jpunpack(  const  unsigned  long  long*  const  Jump_P  )  { 
for  (  int  j=N-l;  j;  j —  )  { 

(  JP  |=  Jump_P[j]  »  32  &  Oxffffffff  )  «=  32; 

(  JP  |=  Jump_P[j]  &  Oxffffffff  )  «=  32; 

} 

(  JP  |  =  Jump_P[0]  »  32  &  Oxffffffff  )  «=  32; 

JP  |=  Jump_P[0]  &  Oxffffffff; 


void  rng::jump(  )  { 

//  compute  jump  state  t.z[0  ..  (N-l)] 
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//  from  jump  polynomial  JP  and  base  states  sy=s.y 

// 

//  t.z[®  ]  =  J[®]*sy[®  ]  +  ..  +  J[d]*sy[d  ]  +  ..  +  J  [Jdmax]  *sy  [Jdmax  ] 

//  t.z[i  ]  =  J[Q]*sy[i  ]  +  ..  +  J[d]*sy[d+i  ]  +  ..  +  J [Jdmax] *sy[Jdmax+i  ] 

//  t.z[N-l]  =  J[Q]*sy[N-l]  +  ..  +  J [d] *sy [d+N-1]  +  ..  +  J [Jdmax] *sy [Jdmax+N-1] 

state  t; 

unsigned  long  long  sy[Jdmax+N];  //  Jdmax+N  states  0  ..  (Jdmax+N-1)  for  jump 
s  =  ref_state;  //  start  s  at  ref  state 

for  (  int  i=®  ;  i< Jdmax+N  ;  i++  )  { 

sy[i]  =  s.y  ;  //  save  states  sy[i  in  0  . .  (Jdmax+N-1)] 

gen()  ; 

} 

memset(  t.z  ,  Q  ,  sizeof(t.z)  ); 
int  i,  k; 

for  (  int  d=©  ;  d<= Jdmax  ;  d++  )  //  JP  term  power  =  d  in  Q  . .  Jdmax 
if  (  JP[d]  ==  1  ) 

//  i  in  Q  . .  (N-l) 

//  k  in  d  . .  (d+N-1)  for  each  d 
//  k  in  0  . .  (Jdmax+N-1)  overall 
for  (  i=©  ,  k=d  ;  i<N  ;  i++  ,  k++  ) 
t.z[i]  "=  sy[k]  ;  //  k  =  d  +  i 

s  =  t  ; 
s.y  =  s.z[Q] ; 
s.n  =  1; 

} 

//  number  of  rngs  instantiated 
unsigned  rng: :ref_counter=Q; 

//  reference  state  for  jump  computation  is  the 
//  zero  state  of  the  previous  instance 
rng:: state  rng : :ref_state; 

//  seed  for  first  rng  instance 
unsigned  long  long  rng : : seed=©xl234Qf®Q ; 

//  Jump  polynomial 

bitsetcrng : : CPdeg>  rng::JP  =  Q; 

//  MT32  jump  coefficients 

unsigned  long  long  rng : : Jump_P[l] [N]  =  {  //  jump  polynomial 

{ 
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//  2"2® 

8x9e48e623e6ef89acull 
8xadfa625443e8e61aull 
8xl3®eb6e39d7ad87cull 
8xbf5e329e®9895a6eull 
8x584d4e62 lb895f 72ull 
8x8ff483831488ald®ull 
8xfl®ff®e88b4581bfull 
8xa®8349e876c2®596ull 
®xcebf®fc45b51bafcull 
8xa6ccfl25®9a68757ull 
®x®318393e711515£6ull 
®xef945eebl6eb2e62ull 
8x7®ffl94®d526ec55ull 
8xf®a4195abf739bl8ull 
8x63£659f499d67dl8ull 
8xf7e392386475cdadull 
8xl932a35ee22974bcull 
8x4 14483 5®a2bl®546ull 
®x6cf®e8a45Qad8d2dull 
®xc58e61afd44aeb3dull 
®xl 149®883cdl846dcull 
8x87229a26d®dbbb61ull 
8x3b98e618ac9463e3ull 
8x91712e41dc®c888eull 
®x627ef®cc8fl52883ull 
8x®d9a8dd2 5  2  77  5bbaul 1 
8x2fc545bl6a57al4cull 
8x5fbbllc5319928beull 
8x9a9a8d54a®662ac9ull 
8xfc7®e7a657deb613ull 
®xb66d®68a69d9de®dull 
®x7953e63abf3®891full 
8x5eb8b7a5dl28cl3®ull 
8x92 1812  If 793cba64ull 
8xle28572dd7a®8888ull 
®x5662faa613f3e®82ull 
®xe4dab®89dc6b857aull 
®xea75dbe3ad33bacbull 
®xe 133 5a5424d33ff lull 
8x57®ae279e6£12517ull 
8x54cb3 lb®8a3®5758ull 
8x9c7618a83887bf8eull 
8x358529c6d991f5cfull 
®xddf46ffl73®dad®6ull 
8xl6®2®3dee6d9511cull 
®x7c6cf7f91bd33c®eull 
®xd456f9beac®d563cull 
®x7dc9dd8699aac3e5ull 
8xl577alc6c47c58a®ull 
®xac6cf9b6f9579eeaull 
®xal®d®fl®dabcf9c9ull 
8x723c7c52c8753946ull 
8xde287f4®e4caba52ull 
®xlell5b4df3bf5befull 
®xb5®5b9eb8bb9ecbfull 
8x741da69a6d96a4c7ull 
Qx458857d6c9cQee6cull 
8xl98bb76e7dl55a28ull 
8x7191cb2138ddede5ull 
8x613d63ca59e57a99ull 
8x48£9c729828acdbeull 
8x2a6bd22a5afce3c8ull 
®xd2 9b 6 7d8 2 4 5 dd4bbul 1 
8xb52d848e88c5fe99ull 
8xa39Q9f3741c39437ull 
8xd61b4b8®67f®4f7®ull 
8x4454acle5325d3f®ull 
8xf952ee4cbb2 3929 lull 
®xb636ff6ba69e4b®bull 
®xea681fafadea8178ull 
8x48f5b4f755d®6e32ull 
®x36®d2fe®75bdbl4full 
8xf7fd81ebe86232c®ull 
®x5a65f6e®bd5da29aull 
8x4c64828a732 IdalSull 
®x22136dba42®d6b6dull 
®xed558e865adf2ea8ull 
8x4112af74a625®®27ull 


8x87bla56b296c3ef9ull 
8xe38e7®67c617c885ull 
8x7f4c4®7b96762a9eull 
8x®532b55cf3a2b23®ull 
8xf536b2d997b48146ull 
®xeefa395aec7c®e7eull 
®xfl®b24af748569b6ull 
®x71a36f46b8f6548eull 
®xf4f9®5544545d811ull 
8x657f4cl7d333c45Qull 
®x®568®88c9625c952ull 
8x357245aabb65e4®8ull 
8x299d5 3792 389 la6aull 
8x8b856886a512b24bull 
8x4b85afl763513442ull 
8x3b3882f8997c4348ull 
8xe®85992912c6be®cull 
8xb7697688b5 187d8bull 
8x697f9b384fl52ca5ull 
8x87d739de4723e6c8ull 
8x5956dac229889£7£ull 
8xe3£a47a938b4b844ull 
8x6918af®39828a7dcull 
8x513edb®bdb83fe58ull 
8xa845 59473 127a429ull 
8x533eea945385f59cull 
8x63fc61abaf68d4b4ull 
8xd95aeb6f®3d2a52cull 
8x®5d®b3538a6f9918ull 
8xa®aaaddbf48c6e68ull 
8xb415a66a8c66aea9ull 
8x7ba41fd®edl9e68bull 
8xf8ee4128529c66alull 
8x679fde4577496188ull 
8xc28b7d76eddda2d®ull 
8x93bedbcd73c29284ull 
8x2 5b289e334846b2 full 
8xb63a741eac495162ull 
8xe4a63af67Q8672d3ull 
8xe87e®4bcb3de®eb®ull 
8xc938192e9f93cd8eull 
8xa®d58658Q32eQ723ull 
8x28bc59474462d4cbull 
8xc52a6dl56bdelb8eull 
8xllcdl426e®a24587ull 
8x3c4c®91f56a8db24ull 
8xa83e868446988782ull 
8xcd89e3eb4fe6ac27ull 
®x79a227af3cc2fc88ull 
8xa488b4696873 laefull 
8xfd8a®36add2b4fbcull 
8xce2b4a8d22bla958ull 
8xe53636528bacal58ull 
®x 1 cbd® 1®  5  5  8  f 3bdadul 1 
®xecdfa3 ld299f 3bf 2ull 
8x2de7a919496bccdaull 
8xb9648ad6el741198ull 
8x2e69a3cl77a94c®eull 
8x9732f889b35b8a91ull 
8xba72d55bla3bc28full 
8xc92988357e3278f3ull 
8x2285fl£c4b2d99a5ull 
8xf9dcl375da31ecaeull 
8xdfa3®45e8®aale®8ull 
8xe9878656cl5£c921ull 
8xae7f9f®5cQ4f9al2ull 
8x7ba91f277da®ba®4ull 
8x8135dl388c82b33eull 
8x8a492e68a94e24f5ull 
8xe56cc7829be38279ull 
8x668392fbb®392cccull 
8xc78e3dl88Q895ac7ull 
8xa34cc64876cl36a7ull 
8xaddf5c2fcc9452c2ull 
8xcaa46®c585d236f6ull 
8x998a57a®b7bbf£feull 
8xd®cd7f485dd7f626ull 
8x®2398e6abfbe9338ull 


8xl®aeda6af61caae®ull 
8x1662 5 3e3322f6158ull 
8xbc8c98818ca42783ull 
8xef®d72875fc84814ull 
8xa®db8ed98ef®b®9cull 
8x59bl369£7e3fdd®8ull 
8x6c553116Q88f3816ull 
8x®eb668b76d58fbf9ull 
8x3aad®ccl9fb3bcc5ull 
8xb76d25a631eaedbeull 
8x8e®a3d31c919626bull 
8x929c958add7e6814ull 
8x2 f24a8abd25ab4f cull 
8xcaa9f68288Q173cbull 
Qxd8f3717e52bd9Q2Qull 
®xf a4edf4ac3 1 177d7ull 
QxdlQ13369728f66a2ull 
8xd5eaec5cc3eblb6cull 
8xae75f9d8d3738734ull 
8xbb62acfe53ee38b5ull 
8x28866ebd8a7a38bcull 
8xdc49472496444e6aull 
8x32a21e®lb32£796bull 
8x52a85dfl389779e9ull 
8x4d8b9e7552d79342ull 
8x2784be952e3549e9ull 
8xd915de®12a6e5b2dull 
8x®®fd99fQfb924d®7ull 
8x973c4d2ce®fad24cull 
8x77429d5283eb965full 
8x2aa41188f2a3®delull 
8x6796637e24f7e®a2ull 
8x6c84f9fe3874e779ull 
8xbd3fe66ae73£297dull 
8xa53bf2647d4befbbull 
8x48bc424d5£79d5eaull 
8x9477 Ifaa8816d93dull 
8xcd48da78558b4a48ull 
8x582a6a393d66b326ull 
8xb362acd897dd8b8aull 
8xe®c5488c2dcebdc4ull 
8xf9c8488f25b5®49cull 
8x5b76e9dd51cc2ad3ull 
8x4c8334e796e27cd4ull 
8xa814dd8dc4bd363dull 
®x3258825®88bleclaull 
8x7df68fb5 122  86c  Hull 
8x29e39e5138353bf6ull 
8x639e7dl3892c943cull 
8xl39a94aade®f47deull 
8x3d®a®fb46al8f587ull 
8xl9ea7£75fdd4f83dull 
8x7a2c656a®d388e28ull 
8x4bf4ea4936acacc7ull 
8xa6ab99174f46e788ull 
8xd7cc3d®34c8c8149ull 
8x919b9fb95Q881152ull 
8xd695d8c6af736e3full 
8x4fe98c4ee®a32c59ull 
®xbae  8ba9557ae231bull 
8xbl3c2cQe84548989ull 
8x7ce8842  58b4Q3Q58ull 
8xa7852757al875822ull 
8xf66c57b57777al9cull 
8x5844c9261966d®elull 
8xl®d39ffed3bb7196ull 
8x585 Ibbdf79e48c8cull 
8xa234a644e8d7715£ull 
8x64a7db5a789e4687ull 
8x776f®99af215d2c®ull 
8xef5cb£47a24b7836ull 
8x64bd58337b538242ull 
8x74311cl448a7b6f8ull 
8xf8®f2281b®cff619ull 
8xf85249abef239557ull 
8xdl73253d4792cb®aull 
8xdfdc2f®73bad8763ull 
8x99d28Q41b®736539ull 


8xal6a4bde28d92e76ull 
8x5856c55495Qb31b£ull 
8x422 1277®d8d49545ull 
8x962288b635®d53e2ull 
8xbb®b®747c388d595ull 
8xbe8b8269219485c2ull 
8xl2d7275a3a5eb5c2ull 
8xl887ebc5cb®c51deull 
8xe4c259df5e33e219ull 
8x88956ca47ad62df9ull 
8xfbe9b9e®b8a69d34ull 
8xlaf3c4b®ld673dcfull 
8xl4b®8ba9ef2458fdull 
8x67feadefc34816f3ull 
8xl5ba53896db7f259ull 
Qx6cQ786efQ8c7aaddull 
8xb457d34fa®88f275ull 
8x6cl7a3decd®ca3eaull 
8x85a24675e7a6aa69ull 
8x382 57895 32e5296cull 
8x2d897c6bfl9b5c5bull 
8xfc76ff3ee66edl23ull 
8xf4fafa74dc2b7884ull 
8xl74684f4b488d9deull 
8x25a327ecldaa9eeaull 
8xe37674ed39blcdc9ull 
8x9cc77871efl9f234ull 
8x8a6bc574a448a366ull 
8xfc39ca9568993d8®ull 
8xd54e695b4771689full 
8xl768cbad7643d4d2ull 
8xfa234881349afb7cull 
8xel9a59ee9f4448a5ull 
8xcal8bb4a6ef67545ull 
8x2ec61ee38blcb88aull 
8xa367ced®3f8f®278ull 
8xb4f76bd83a673955ull 
8x38e84f58f4cadd2cull 
8x85cbddb22eddd4d3ull 
8x547d®63de839f589ull 
8x7el238facc7356d7ull 
8xble®lda354e632alull 
8x5b7£dl381b766c23ull 
8x6ab3b24874552®29ull 
Qxe®99e812Q18b88e4ull 
8x7e833aacdccl455cull 
8x5bdfd4297a37965aull 
Qx238fl67dc8ea®6Qbull 
8x787e77f8ac8ee42bull 
8x8 14e62 Ide9e6c2b9ull 
8xfd3413e3c3a32bcbull 
8x6cfd9623b74ecf92ull 
Qxcd43e9Q4a2Q26126ull 
8x892flcfe318192dbull 
8xe61d®15825118648ull 
8xe758a®b4a658®f2eull 
8xl38d58687b84db2full 
8x6bcd3d9d81eaabc7ull 
8x4c58d®26a®e38e27ull 
8xb2cc7386cf967d26ull 
8x65ebcdd5392d4d46ull 
8x4559fae2c9573f®lull 
8xf4a®6ae488fa®4faull 
8xe6348bf7dcc7cdc4ull 
8x9654291782475d4cull 
8xc4a37afl3746625aull 
8x725e3d79b®ce59f2ull 
8xcblb2b®848fce64eull 
8xf597e6d649379b9£ull 
8xd®6c8ec8ad7dc968ull 
8x4ce3dc55f68ee4b5ull 
8xf75fe8f3c£lceb®£ull 
8x99487bl7a582d52eull 
8xe6a58d7189e7e77®ull 
8x639ffded4c8c®b6£ull 
8xe386f®18ea23c96cull 
8x535372dddlcd®dl5ull 
8x8Q®Q8Q8Q44337dd9ull 
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//  2' 1SSQ® 

//  Qxe4115b36993bffe7ull 
//  Qxfff75a273cd88961ull 
//  Qxf7fc813dd4bae594ull 
//  Qxea3b8b69d445696bull 
//  Qxf54d97e®ebel8662ull 
//  Ox6cb3b4ee2f7e®de5ull 
//  ®xbel9ae26ca4£53®8ull 
//  Qxl®f7®a76dbc6ad2eull 
//  8x9fd856495627c455ull 
//  8x85115a2a4452a8f6ull 
//  ®xd5331a®85®bf898®ull 
//  8x515a7Q2d83b7658dull 
//  8x73dde418b88e53a5ull 
//  ®xd2  7 f 2  6ada47£3  7f ®ul 1 
//  ®x®9e 1 7 648d®7®6b®aull 
//  8x52fd314c5®a6b3e7ull 
//  8x867 19cebc®b84b6cull 
//  8xa®5b8 17285 5®7b5cull 
//  8x®e4a7dl635e6elldull 
//  8xff2a®f842®af8ef4ull 
//  8xe®b57515e8bc958aull 
//  8x25cd43bc7b57c283ull 
//  8xd83b®34ed5®84c98ull 
//  8x47bc91a®515eb478ull 
//  8xd861af®c®lbdb29eull 
//  8xlbf4733e638871baull 
//  8xe67el7c3338c92c8ull 
//  8x878363 5 5bd4d79e5ull 
//  8x792f81f225e915dcull 
//  8xfbc®7978eb393845ull 
//  Qx®84b91abl8c2a9efull 
//  8x768282 3d2a49464eull 
//  8x44db4e®532d8518bull 
//  8x58f821f4615a7£d5ull 
//  8x998535d74126aebfull 
//  Qxc8be33a8fffd®2f2ull 
//  8xff32fc®bc88c36b5ull 
//  8x43b8bc88be3ele®eull 
//  8x3cf391ac21ceded5ull 
//  8x94f9c8bbd82bl2feull 
//  8x382 39324f®66d2baull 
//  8x3748faf2e®ecaa35ull 
//  8xa25522abbb4e6583ull 
//  8x48d89ab57e®6Q8c£ull 
//  8xbda6be36455bfb62ull 
//  8xle595b93e437a56eull 
//  8xb65£114ae2c72774ull 
//  8xllla3582ac2391alull 
//  8xb7c4863879db7cf7ull 
//  8x8adc73fd6dl7balbull 
//  8x3®ddf2fefalflbeeull 
//  8xfb987e989718fdd3ull 
//  8x559d58a®c9cc35dlull 
//  8x6659fe685®67eab4ull 
//  8x35f2ab22d3e®7eb£ull 
//  8x89a®4e58ffaal242ull 
//  8x8138235729aaea77ull 
//  Qx99Qc 1 lb 1 9bd74bcaull 
//  8x4772ebd589571283ull 
//  8x35af®ldl793eabb6ull 
//  8xfe784e81a58b4684ull 
//  8xel53193e61d62943ull 
//  8xdf 5  7ba®bQ2  2  7 8bb®ul 1 
//  8xe4446755d®acd36aull 
//  8xb6dfb3226f22d513ull 
//  8x31fc6248e3182d8cull 
//  8x8724aal2a5b781deull 
//  8x3a54e88®d7e42e8bull 
//  8x62533fc®7122885dull 
//  8x6elae218b67a4499ull 
//  Qxfb234c25afccQ8b7ull 
//  8x624856becl6ec681ull 
//  8xlf3d9de67c3f8354ull 
//  8x6b79ee94dd67be2bull 
//  8xcebll®b52b78449bull 
//  Qxf®5895253177fec4ull 
//  8x272 Ifdb614a2831full 
//  8x26b£f7143c5d688eull 


8x2298d948e83a43acull 
8x676783d88a3bf981ull 
8x544d®855917c2fdlull 
8x391aeaffcad482e3ull 
8x63869eee5f596911ull 
8xf389adbdd422®ccdull 
8x®79baed9cd®986a4ull 
8x5a5 lb2 12b9a®83 lbull 
8x8cecae753e7428adull 
8xfc9bc3 16f 3b®2 5e4ull 
8xl93ce2885edc6c®eull 
8xc66bc86dl4fld331ull 
8x®52ba71c4c726681ull 
8x8893e3£5eafcb587ull 
8x4e27186995e7cd98ull 
8xe8456586297599edull 
8xb5f®693a9d64652bull 
Qx5ef®55acd7fa®c82ull 
8xf9ab85ae5dalc®c6ull 
Qx7ddd88bla®Q3ccl4ull 
8x4eb2ebla®793912eull 
8x98b7e7224988918bull 
8x5b641c£2fbf9ebf7ull 
8xcb878cc53fad6c93ull 
8x6c8c898b5fdc9fccull 
8xb7811ffb6f4435c9ull 
8xcla®ee38795c6d75ull 
8x38d467b9274591e6ull 
8xl2ef3163247d7379ull 
8x6dl441a3f9851dc3ull 
8x557876d5e4d9ac7dull 
8xdfb®b41fd3982e38ull 
8x717b47fb£78d5c5full 
8x37486ca8e8d46bcdull 
8x63d®7cl7£a2a4258ull 
8x3e293178ea7ffc8dull 
Sxaae  5baabb8af 484  5ull 
Qx41ef®8cfff38Q9Qeull 
8xe45dfcdc5982®adaull 
8xf8dd6e®e3e82cda6ull 
8x489b5832e4ae763aull 
8x3585a46928686b95ull 
8x68e583ddffdabe4eull 
8x3453 5 l£c2138979cull 
Qx2c8£875cc9£5 1 5®2ull 
8x2 1 5bba®8bea® ld8eull 
8x919889fd8blace2aull 
8xdl56382ef®71£675ull 
8xab38125457df2e67ull 
8x27ac93aaa8ab83feull 
8x91eff®deel4f8c6cull 
Qxea33 1597241 122 5eull 
8xblb4bcc4cl®12f5dull 
8x69882 3774 lfc96beull 
8x38d®f32872a8eba4ull 
8xa®987c638e73ec62ull 
8x26361a872163812cull 
8x2fb2b6a748cc74d7ull 
8x89894e3253332726ull 
8xce83 17862 17b7abcull 
8xdb7215al£d£38984ull 
8x4e986c7b876264a8ull 
8x84bf35dlld58896cull 
8xd23f9baaff®3®2a®ull 
8xl£ce3748689a4889ull 
8xecfc9d35c95a953dull 
8x4dd48855cb8398f7ull 
8x5d8d44e92f3fbec8ull 
8x235b72618314d2e6ull 
8x2dd68c417df6e2a6ull 
8xffc6b5c3641aff53ull 
8x787all775bf®645eull 
8xl87adaae34d97blcull 
8xe8e37f 77£5 19faceull 
8xf7f488c5911c538dull 
8x7£d533f6bd57f96cull 
8x2e®da59be7433f6bull 
8xa3cff887el57c93eull 


8x2277el2bc9c®74b4ull 
8xle833b82c93e2fQbull 
Qx7c46429c52e®d4adull 
8x37fc6®7c3c5f®e96ull 
8xb8e381dl782841abull 
8x7e628e6e766®b595ull 
8xd2736da73c95794full 
8x7cl4e6393859c®3aull 
Qxl 5 5493428 17d2da®ull 
8x873 5faa®187a7eadull 
8xc8927fd78ffl3d62ull 
8x918bea7b497abaf4ull 
8x®95ee33eddc88352ull 
8x2b37a8de97125d®bull 
8x6e26f66817ef329dull 
8x22al943a£e2e2761ull 
8xfb34e899be436b5aull 
8x3c66cd88c389d684ull 
8x835bbebf54cf56ffull 
8x935ed3846eb®69e®ull 
8x2cbd8ce®5a2cc449ull 
Qx719cab883b5cl7d3ull 
8x4b588al95cl38beaull 
8x9a6474c53da37aa3ull 
8xcle9e98123dccdb8ull 
8x4  fade  7c 1 5  5®bdf 2  Sul 1 
8xlb41Qe29888de987ull 
8x711fe3bc4a56c662ull 
8xbd5fb88bad®ba338ull 
Qx22d599blda7b9ee5ull 
8x4831£c382261acaeull 
8x67f8325a6a38152cull 
8xd9d452fel9698c77ull 
Qxc4a9ffc7e3£525£5ull 
8xf348876ea®d48ba5ull 
8xe3cfef2f7cee9c98ull 
8x5882b7cc6albaaa9ull 
8xd5819f4b24a896a®ull 
8xee442aeb69b6bl85ull 
8xcee3e®e6fda76814ull 
8x514bddab718ee7ecull 
8xlc8eb2 168a63c679ull 
8xf41f2ac3£1222373ull 
8x37264 1£56® If 65f6ull 
8xdaa82  5db®dad®7alull 
8x£5578592d31532ccull 
8x81cb2b8de®c2e998ull 
8x3714ba5d26al75e3ull 
8x7fd44d82c8ca2 145ull 
8xa5c®988af55ea944ull 
8x®b5bbf68f6bf8e8dull 
Qxa®f26aba7293d313ull 
8xf6a64a4b668159d6ull 
8x®237386c4c£e7ac7ull 
8x77c276bae36b9fe6ull 
8x4181d4491827e5baull 
8x85595a36b®6dallfull 
8xe8592fdcc4cfd658ull 
8x563883733 3£3fb98ull 
8xd76192e68bQdb9Qdull 
8xl9ffed6dalddl6elull 
8x569blb2e61c2bb66ull 
8x39cab41e461ee546ull 
8x85294873c®51cd53ull 
Qx9ad®a97eeQ611a6full 
8xcc51c®865c76dc3aull 
8x®e42b7£914de69bdull 
8xc888bb32aabce373ull 
8xe2b74c331b4de®deull 
8x3®c®888d®ee4®6baull 
8x8e5b31ef54Q55d29ull 
8x39bf99598959ed5bull 
8xb8f 5dd4ec2  34a3  5®ull 
8x28539659c39253f®ull 
8x92afae9e8c7c81fcull 
8xc58fclbd74828485ull 
8x2a®99cc7588a726aull 
8xbe88f3d6acb68f48ull 


8x33cl9b98elb9deadull 
8x4f®9b86e86159fd9ull 
8x7b82ef96b5e8b98aull 
8x4e6ff83ca4ffe513ull 
8xad664484a638f®ldull 
8x659deal5686f6aaeull 
8x91eab4989d8flce2ull 
8xe7ff466589194169ull 
8x61e3fcal66da2712ull 
8xe48d5396b2ef2a98ull 
8xl92826468ae8b®c5ull 
8x2 13b445 12dl39al7ull 
8x45bdb5ea5a88a229ull 
8xb®f2 139f 78dd58ddull 
8x4aa4ad7e442865b2ull 
8x77e3827®18b859b7ull 
8x3297798786df®b35ull 
8xe889f 6c5e®cf2  5  5  5ull 
8x43763a856ea®a9b4ull 
8xlae®6586dd33c79full 
8x48f7bbdbdb®4ed8aull 
8xd267ddd5d®fle45full 
8xf®e775a546bb9b8cull 
8xl68b8f4418c77948ull 
8xfe®1889db698de64ull 
8x239decb2d26fede2ull 
8x8c62bab3b927dd5full 
8x727959d6f65b3alfull 
8xb7e7f®436188c4dbull 
8xef®b48f®8clb9438ull 
8xb5dad81bd234b39eull 
8x6113d49dba637bdeull 
8x4d2849a4c44a4c81ull 
8xd8a®c62f35£la934ull 
8x938845349d782c3dull 
8xlc3ad®8®5a53d5caull 
8x2 17d4dc96a4 If 3 1 lull 
8xd8ee®ldfaf4845c6ull 
8x514aab282a®414alull 
8xee®ccf686463c893ull 
8x227385bdbl®7c84bull 
8x9f624dabl5c58a7®ull 
8xac88782bddcafc2dull 
8xa746c8dl®b5 Ie6e5ull 
8x6dc7e7cb53616925ull 
8x3beb9b5bla®bda28ull 
8xf7891c58a42c4229ull 
8xf8f®4fad®c® 1489 lull 
8x25147538a78ccc4full 
8x41f5c®8cdb6997e6ull 
8x533ba3687e45cbc7ull 
8xfbafb565298dc758ull 
8x8295deea6cb8dc35ull 
8x94784af94e®51clcull 
8xbb326353fcba584cull 
8xb3d48bc®7c91ce8full 
8x6d5cd34a95fbl355ull 
8xa3272®74®82 lde5dull 
8x818b387dc21da55eull 
8xdbec797el4385e71ull 
8xe979bl®a385e75blull 
8x72f981e988fa®dc3ull 
8x9a5a8f69adla2b3bull 
8xlc687c62af6767b9ull 
8xab438f 8e682  5b4baull 
8xae8582bddcb487abull 
8x3993881ad83feQ67ull 
8xcc98dfcafa49f98full 
8x®b41b7abba3a4cdaull 
8x824d6Qd26c8d3375ull 
8xd4197118e5e7®bb3ull 
8xebab888bb8be®182ull 
8xda6 le  5  4cdcdf 2  2  f 9ul 1 
8x8f23f25a®eab6773ull 
8x77cc28e3468db298ull 
8xd71e77®f®8ae®746ull 
8x3Qd38b3695®c3®8Qull 
8x888888® Ib3fdecd4ull 
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NO.  OF 

COPIES  ORGANIZATION 


1  DEFENSE  TECHNICAL 
(PDF  INFORMATION  CTR 
only)  DTIC  OCA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FORT  BELVOIR  VA  22060-6218 

1  US  ARMY  RSRCH  DEV  & 
ENGRG  CMD 
SYSTEMS  OF  SYSTEMS 
INTEGRATION 
AMSRD  SS  T 
6000  6TH  ST  STE  100 
FORT  BELVOIR  VA  22060-5608 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
IMNE  ALC  IMS 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  T 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 


ABERDEEN  PROVING  GROUND 
1  DIR  USARL 

AMSRD  ARL  Cl  OK  TP  (BLDG  4600) 
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NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


AMSRD  ARL  SL  BA 
D  FARENWALD 
M  PERRY 

AMSRD  ARL  SL  BD 
R  GROTE 

AMSRD  ARL  SL  BG 
S  SNEAD 

AMSRD  ARL  SL  BS 
W  WINNER  (4  CPS) 

AMSRD  ARL  SL  BW 
L  ROACH 


ATRC  W 
A  KEINTZ 

WSMRNM  88002-5502 

1  USARL 

AMSRD  ARL  SL  E 
R  FLORES 

WSMR  NM  88002-5513 
ABERDEEN  PROVING  GROUND 

1  US  ARMY  DEV  TEST  COM 
CSTE  DTC  TT  T 
3 14  LONGS  CORNER  RD 
APG  MD  21005-5055 

1  US  ARMY  EVALUATION  CTR 
CSTE  AEC  SVE 
R  LAUGHMAN 
4120  SUSQUEHANNA  AVE 
APG  MD  21005-3013 

15  DIR  USARL 

AMSRD  ARL  SL 
J  BEILFUSS 
J  FEENEY 
J  FRANZ 
M  STARKS 
P  TANENBAUM 
AMSRD  ARL  SL  B 
G  MANNIX 


1  ASST  SECY  ARMY 
(CD  ACQSTN  LOGISTICS  &  TECH 
only)  SAALZPRM2E661 

103  ARMY  PENTAGON 
WASHINGTON  DC  20310-0103 

1  ASST  SECY  ARMY 
(CD  ACQSTN  LOGISTICS  &  TECH 
only)  SAAL  ZS  RM  3E448 

103  ARMY  PENTAGON 
WASHINGTON  DC  20310-0103 

1  DIRECTOR  FORCE  DEV 
(CD  DAPR  FDZ 
only)  RM  3A522 

460  ARMY  PENTAGON 
WASHINGTON  DC  20310-0460 

1  US  ARMY  TRADOC  ANL  CTR 
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NO.  OF 

COPIES  ORGANIZATION 


ABERDEEN  PROVING  GROUND 


21  DIRUSARL 

AMSRD  ARL  SL  BA 
G  BRADLEY 
R  DIBELKA 
M  KUNKEL 
AMSRD  ARL  SL  BD 
J  COLLINS  (5  CPS) 
A  DRYSDALE 
T  HOLDREN 
L  MOSS 
E  SNYDER 
AMSRD  ARL  SL  BE 
K  BATES 
R  SAUCIER 
AMSRD  SL  BS 
J  ANDERSON 
J  AUTEN 
R  BOWERS 
M  BURDESHAW 
E  DAVISSON 
V  PHU 

G  SAUERBORN 
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Intentionally  left  blank. 
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